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SECTION  I 


INTRODUCTION 


1.0  BACKGROUND 
1.1  Orientation 

Because  of  their  Inherent  flexibility  of  control,  reliability,  and 
efficiency,  uses  of  electromechanical  actuators  and  electric  drives  are 
finding  increased  interests  in  aircraft  applications.  Three  areas  of 
technology  advancement  over  the  past  decade  are  largely  responsible  for 
placing  electromechanical  energy  converters  into  a  favorable  position 
when  compared  to  hydraulic  motors  or  actuators: 

(a) .  Development  of  rare  earth  permanent  magnet  (PM)  motors  with 

Inherent  high  efficiency  and  high  power  to  weight  ratios. 

(b)  .  Advancements  in  power  level  solid-state  devices  with  high 

switching  speeds. 

(c) .  Emergence  of  microprocessors  allowing  control  capability 

and  sophistication  far  surpassing  that  of  analog  systems 
while  reducing  the  volume  occupied. 

The  conventional  or  commutator  DC  machine  performance  characteristics  in 
the  areas  of  speed  control  and  position  control  are  highly  desirable. 

The  brushless  DC  machine  has  wide  range  torque-speed  characteristics  like 
unto  the  commutator  DC  machine  without  the  commutator-brush  maintenance 
problems.  In  addition,  the  brushless  DC  machine  with  a  permanent  magnet 
rotor  has  certain  superior  features  to  the  conventional  DC  machine: 

(a) .  Field  excitation  is  eliminated,  which  removes  the  complexity 

of  supplying  power  to  a  rotating  member.  Also,  machine 
efficiency  is  increased  due  to  absence  of  field  excitation 
losses . 

(b) .  Higher  speed  design  is  possible  for  PM  rotors  than  is  feasible 

with  wound  rotors  permitting  increased  gear  ratios,  which  leads 
to  substantial  reduction  in  electric  machine  power  to  weight 
ratios . 

(c)  .  Thermal  transfer  characteristics  are  improved  .since  the  bulk 

of  the  losses  (ohmic  and  core  losses)  are  generated  within  the 
stationary  member,  allowing  efficient  implementation  of  fluid 
cooling. 


1.2  Identification  of  Problem 

Brushless  DC  motor  performance  is  reported  in  the  literature,  but  it 
concentrates  on  the  nature  of  instantaneous  voltage  and  current  along  with 
average  values  of  developed  torque  Cl -9] .  These  works  typically  describe 
systems  that  do  not  operate  continuously  at  near  zero  speed,  and  thus, 
only  average  value  of  torque  is  of  concern.  However,  the  brushless  DC 
machine  inherently  has  oscillatory  components  of  instantaneous  torque 
at  all  frequencies  that  are  integer  multiples  of  six  times  the  electrical 
radian  of  the  stator  impressed  voltage.  Since  stator  frequency  is  directly 
related  to  rotor  position,  these  oscillatory  components  are  in  the  fre¬ 
quency  range  of  mechanical  system  response  at  low  speed  values. 

Published  works  relating  to  performance  of  brushless  DC  machines 
have  not  analyzed  cases  of  position  or  sustained  near  zero  speed  operation; 
thus,  pulsating  torque  components  have  been  considered  of  no  significant 
consequence.  However,  Williamson  et  al  ClOD  have  mentioned  that  torque 
capability  deteriorates  at  low  speeds,  and  Demerdash  and  Nehl  have  shown 
some  instantaneous  developed  torque  and  power  wave  forms  Cll3  without 
comment  on  the  oscillatory  components. 

Widespread  usages  of  brushless  DC  motor  drive  systems  in  low-speed 
and  position-control  actuator  applications  are  contingent  upon  development 
of  control  philosophies  and  hardware  realizations  of  power  conditioning 
arrangements  that  allow  bidirectional  power  flow  while  resulting  in 
Instantaneous  motor  developed  torques  that  are  free  of  harmonics  in  the 
range  of  response  for  coupled  mechanical  loads. 

1.3  Consideration  of  Alternative  Technologies 

There  are  two  basic  power  electronic  circuits  that  are  used  as  power 
conditioning  links  to  couple  brushless  dc  motor  drives  to  a  high  frequency, 
three-phase  AC  source: 

(a)  DC  link  inverter 

(b)  Cycloconverter  link 

Power  circuits  of  these  two  basic  approaches  are  illustrated  in  functional 
block  form  by  Figure  1.  Variations  of  each  arrangement  are  made  .according 
to  the  control  philosophy  adopted  to  satisfy  required  motor  performance. 
Further,  for  low  voltage  (500  V  or  less),  low  current  (200  A  or  less) 
applications,  power  circuits  can  be  synthesized  with  transistors  for 
controlled  switching  elements,  while  for  high  power  level  applications 
silicon  controlled  rectifiers  (SCRs)  must  be  used  for  controlled  switching 
elements.  For  this  study,  SCR  switching  elements  are  assumed. 

1.3.1  Operation  for  Speeds  Not  Hear  Zero.  For  the  case  of  a  DC  link 
inverter  driving  a  brushless  DC  motor  that  does  not  operate  at  near  zero 
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Figure  1.  Alternative  Power  Circuits 

(a)  DC  Link  Inverter 

(b)  Cycloconverter  Link 


speed,  the  rectifier  may  be  a  diode  bridge  while  the  chopper  is  used  to 
vary  magnitude  of  DC  voltage  applied  to  the  inverter  input  terminals;  or, 
the  chopper  may  be  eliminated  and  a  phase-controlled  converter  used  as 
the  rectifier  to  vary  inverter  input  voltage. 

For  the  case  of  a  cycloconverter,  the  necessity  of  a  DC  link  does 
not  exist,  and  thus,  there  is  only  one  power  conditioning  module  inter¬ 
connecting  the  three-phase  AC  source  with  the  brushless  DC  motor. 

Generally,  the  cycloconverter  may  be  realized  as  either  a  midpoint  or  a 
full-bridge  arrangement  and  either  phase-control  or  synchronous  envelope 
operation  implemented. 

1.3.2  Operation  for  Speeds  Near  Zero.  For  near  zero  speed  operation  of 
a  DC  link  inverter  for  high  power  level  systems,  the  counter  emf  of  each 
motor  phase  is  small  enough  so  that  natural  commutation  of  the  inverter 
SCRs  is  prohibited.  In  such  a  case,  the  chopper  is  used  to  reduce 
inverter  terminal  voltage  to  zero  sufficiently  long  to  accomplish  commu¬ 
tation  of  the  inverter  SCRs  C3] .  Since  the  chopper  must  be  controlled 
so  as  to  enhance  inverter  commutation,  the  rectifier  must  be  a  phase- 
controlled  converter  to  permit  necessary  inverter  input  voltage  magnitude 
control. 

For  near  zero  speed  operation  of  a  phase-controlled  cycloconverter 
using  circulating  current  free  mode,  discontinuous  load  current  tends  to 
occur,  leading  to  Increase  in  load  current  harmonics.  Continuous 
current  can  be  restored  by  changing  to  circulating  current  mode  of  control 
at  the  expense  of  increased  losses  due  to  circulating  reactive  current  C123. 
Control  is  generally  simplified  at  fractional  hertz  frequency  operation  by 
use  of  synchronous  envelope  control  with  circulating  current  free  mode. 

1.3.3  Operation  for  Regenerative  Power  Flow.  For  the  case  of  a  DC  link 
inverter,  regenerative  power  flow  requires  further  power  circuit  modifica¬ 
tion.  Gating  signals  to  inverter  SCRs  are  suppressed  and  the  associated 
shunting  diodes  function  to  form  a  full-bridge  diode  rectifier.  A 
switching  circuit  is  introduced  in  the  DC  link  to  reverse  polarity  of  the 
voltage  appearing  at  the  DC  terminals  of  the  rectifier  module.  The 
rectifier  must  be  a  phase-controlled  converter  operated  in  the  synchronous 
inversion  mode.  Reversal  of  motor  developed  torque  requires  three 
coordinated  control  actions:  phase  forward  and  suppression  of  inverter 
SCRs;  polarity  reversal  switch  activation;  and,  phase  forward  of  phase- 
controlled  converter  SCRs. 

Regenerative  power  flow  using  the  cycloconverter  link  is  introduced 
by  simply  delaying  SCR  firing  angles  beyond  90°.  Rate  of  change  of  motor 
developed  torque  from  a  positive  to  a  negative  value  is  controlled  by  the 
rate  at  which  the  SCR  firing  angles  are  changed.  Thus,  the  nature  of 
transition  from  motoring  to  regeneration  is  determined  by  a  single  control 
action  leading  to  smooth  change  with  minimum  time  delay. 


2.0  OBJECTIVE 


Of  the  presently  available  power  conditioning  technologies  for 
linking  the  brushless  DC  motor  to  a  high  frequency  polyphase  AC  source  in 
a  near  zero  speed  or  actuator  application,  the  cycloconverter  using  syn¬ 
chronous  envelope  control  with  circulating  free  current  mode  of  operation 
appears  most  promising  in  that  it  requires  no  external  commutation  cir¬ 
cuitry  and  it  can  be  smoothly  and  quickly  changed  from  motoring  to 
regeneration  by  a  single  control  action.  This  arrangement  using  a  mid¬ 
point,  three-pulse  cycloconverter  is  chosen  for  study  with  pulse  width 
modulation  (PWM)  control  to  eliminate  those  harmonics  from  the  output 
waveform  that  lead  to  harmonics  in  the  developed  torque  in  the  range  to 
which  coupled  mechanical  loads  can  respond. 

3.0  SCOPE 

The  following  tasks  were  carried  out  to  accomplish  the  objective: 

(a)  .  Determine  an  appropriate  mathematical  model  of  the  system  for 

near  zero  speed  operation. 

(b)  .  Determine  the  current  harmonics  to  be  eliminated  in  order  to 

remove  undesireable  pulsations  in  motor  developed  torque. 

(c) .  Develop  a  modulation  function  to  use  in  control  of  phase 

voltage  for  eliminating  the  undesired  current  harmonics. 

(d)  .  Calculate  motor  performance  over  a  wide  speed  range  with  and 

without  the  PWM  phase  voltage. 

(e)  .  Compare  ohmic  losses  with  and  without  PWM  to  assess  the  effect 

of  harmonic  elimination  on  efficiency. 
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SECTION  II 


PROCEDURE 


SYSTEM  DESCRIPTION 


Typical  values  for  motor  parameters  were  obtained  by  ratio  of  known 
values  reported  In  the  literature  Cll3  where  the  results  are  given  In 
Table  1.  The  motor  is  wye  connected  with  a  magnum  design  speed  of 
45000  rpm. 


TABLE  1.  MOTOR  PARAMETERS 


Parameter 


Symbol 


Value 


No.  poles 
EMF  constant 
Stator  resistance 
Stator  leakage  Inductance 


0.0225  V- s/rad 
0.4  0 
25  pH 


Characteristics  of  the  balanced  three-phase  source  are  listed  in 
Table  2. 


TABLE  2.  SOURCE  CHARACTERISTICS 


Quantity 


Value 


Phases 

Line  voltage 
Frequency 


3 

138  V 
7950  Hz 


The  power  circuit  schematic  of  a  midpoint,  three-pulse  cycloconverter 
linking  the  high  frequency,  three-phase  AC  source  to  the  brushless  DC 
motor  is  depicted  by  Figure  2.  Use  of  the  switch  in  the  midpoint  or 
neutral  line  will  be  discussed  later. 


ASSUMPTIONS 


When  analyzing  rare  earth  permanent  magnet  machines  with  non-magnetic 
retaining  rings,  it  has  been  found  that  position  dependence  and  interphase 
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coupling  of  machine  inductances  can  be  neglected  CllD.  Such  an  approxima¬ 
tion  leads  to  decoupled  equations  that  can  be  used  in  networks  formed  by 
addition  of  power  conditioning  circuitry  with  minimum  difficulty.  Motor 
counter  emfs  are  taken  to  be  sinusoidal.  Further,  since  the  study  of 
this  report  is  concerned  with  low  speed  operation,  the  commutation  overlap 
is  very  small  compared  with  a  period  of  the  motor  current  wave  form  and 
can  be  neglected. 

The  SCRs  are  modelled  by  a  0.02  ohm  forward  resistance  (2  V  at  100  A) 
and  a  blocking  resistance  described  by  1000  i2  ohms. 

3.0  MATHEMATICAL  MODEL 

Referring  to  Figure  2,  using  the  above  assumptions,  and  applying 
Kirchhoff's  voltage  law  results  in  the  following  set  of  simultaneous 
differential  equations  to  describe  the  motor  electrical  performance: 

r  ('11R  +  van  *  el>  '  <« 

La 

~  <-i2*  +  vbn  -  e2)  (2) 

La 

^  <-*3*  +  *cn  -  «3>  «) 


Ri,  R2,  and  R3  are  the  sum  of  the  motor  phase  resistance  and  the  resistance 
of  the  particular  SCR  between  the  source  and  phase  of  the  motor  at  the 
particular  Instant  of  time  under  analysis.  The  form  of  phase  voltages 
(van»  vbn»  vcn^  be  discussed  later. 

Neutral  connection  current  is  given  by  Kirchhoff's  current  law  as 


*nN  "  +  *2  +  <*> 

As  long  as  motor  speed  (w0)  remains  non  zero,  the  electromechanical 
developed  torque  is  given  by 


Td  “  (eiix  .+  *2*2  +  e313)/wm  (5) 

Equations  (1)  -  (5)  completely  describe  the  instantaneous  performance 
of  the  brushless  DC  motor. 


di 


1 

dt 


" 2 
dt 


di3 

“dt 


4.0  CURRENT  HARMONICS  TO  BE  ELIMINATED 


For  both  the  case  of  a  DC  link  inverter  drive  and  the  case  of  a 
synchronous  envelope  cycloconverter  drive,  the  motor  phase  currents 
approach  the  pulse  width  controlled  square  wave  of  Figure  3.  A  Fourier 
series  representation  of  the  wave  form  shows  that  all  odd,  nontriplen 
harmonics  exist  in  the  phase  current. 

i.  *  —  T  —  cos  (— )  sin  nut,  n  *  1,5,7,11,13,  ••  •  (6) 

1  it  Z.  n  6 


Figure  3.  Approximate  Form  of  Motor  Phase  Current 


Similarly,  the  remaining  balanced  phase  currents  are  given  by 


i3  ■  ^  I  j  cos  (HI)  Sin  n(ut  +  )  (8) 

v  “  n  6  3 

n  •  1,5,7,11,13, • • • 


Since  the  motor  counter  emfs  (e^ej.e,)  are  sinusoids  of  fundamental 
frequency  and  form  a  balanced  three-phase  set,  each  can  be  expressed  as 


10 


Ejjj  sin  (wt  +  $) 


(9) 


e2  “  ®m  8in 


e3  -  ^  sin  (wt  +  <j>  + 


(10) 

(11) 


If  equations  (6)  -  (11)  are  substituted  into  (5),  the  simplified 
result  for  instantaneous  developed  motor  torque  is  found. 


4  Enl  l  i  cos  (HI)  {cosC(n-l)I  wt  -  $3  -  cost(n+l)I  wt  +  4*3>(12) 
"  n  n  o  o  6 


From  (12) ,  it  is  seen  that  the  instantaneous  developed  motor  torque 
is  made  up  of  a  constant  term  plus  all  multiples  of  the  sixth  harmonic 
of  motor  phase  current  frequency.  Further,  the  multiples  of  sixth 
harmonic  components  of  torque  are  a  direct  result  of  the  current  harmonics 
that  lie  immediately  on  either  side  of  multiples  of  six.  Specifically, 
the  sixth  harmonic  of  torque  results  from  the  fifth  and  seventh  harmonic 
of  current;  and,  the  twelfth  harmonic  of  torque  results  from  the  eleventh 
and  thirteenth  harmonic  of  current;  etc.  Clearly,  the  motor  developed 
torque  will  be  constant  if  the  harmonics  of  motor  phase  current  described 
by  6n  ±  1,  n  -  1,2,3,-..  are  eliminated. 


The  above  discussed  harmonics  only  need  to  be  eliminated  if  their 
existence  results  in  torque  harmonics  in  the  range  of  response  of  a 
coupled  mechanical  load.  Since  typical  mechanical  loads  exhibit  negligible 
response  to  frequencies  above  20  hertz,  it  should  only  be  necessary  to 
eliminate  those  nontrlplen,  odd  current  harmonics  beyond  the  first  for 
which 


M 


80  it 
6  P“m 


(13) 


where  p  ■  number  of  poles  and  wD  is  the  motor  speed  in  rad/s. 


5.0  MODULATION  OF  PHASE  VOLTAGE 


The  cycloconverter  system  is  voltage  excited;  thus,  if  the  (n  ±  1) 
harmonics  are  to  be  eliminated  from  motor  phase  currents,  they  must  be 
eliminated  from  the  phase  voltage  set. 

5.1  Development  of  Modulation  Function 

Modulation  function  h(wt)  that  contains  no  selected  harmonics  up 
through  M  is  derived  in  Appendix  B.  The  balanced  three-phase  voltage  set 
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to  be  applied  for  equations  (1)  -  (3)  is  given  by 


van  " 

vd  h(ut) 

(14) 

vbn  “ 

2tt 

vd  h(ut - j  ) 

(15) 

vcn 

Vd  h(o»t  +  —  ) 

(16) 

If  is  Che  radian  frequency  of  Che  three-phase  AC  source,  Chen  vd  is  a 
periodic  wave  form  described  by 


a 

s 

ft 

-  Vffl  sin(a)Lt  +  1  +  ag)  ,  0  <  a^t  < 

(17) 

and,  vd(t  +  T) 

"  vd(t) 

(18) 

where  T 

2tt 

3ml 

It  is  shown  in  Appendix  A  that  if  Che  modulated  phase  voltages 
(van»vbn»vcn)  a*e  to  contain  none  of  the  selected  harmonics  through  M, 
it  is  necessary  that  the  frequency  of  vd  be  much  greater  than  the 
frequency  of  h(t);  however,  this  condition  is  easily  satisfied  at  near 
zero  speed. 

Figure  4  illustrates  the  unmodulated  phase  voltage.  The  modulated 
phase  voltage  is  shown  by  Figure  5. 

5.2  Minimization  of  Slack  Variable 

Development  of  the  modulation  function  for  use  with  the  synchronous 
envelope  cycloconverter  in  Appendix  A  required  an  extension  of  the  well- 
known  pulse  width  modulation  techniques  for  elimination  of  harmonics  in 
inverter  drives  Cl4]  to  include  a  slack  variable.  The  modulation  function 
contained  no  selected  harmonics  through  M  as  determined  by  (13) .  For 
the  reported  work,  the  slack  variable  is  always  chosen  so  as  to  minimize 
the  sum  of  the  squares  of  the  amplitude  of  the  next  higher  (6n  ±  1) 
current  harmonic  pair  beyond  harmonics  selectively  eliminated.  Listings 
of  the  FORTRAN  programs  used  to  compute  the  modulation  function  are 
presented  in  Appendix  B. 


5.3  Cycloconverter  Connection  Method 

Since  the  modulation  function  h(uit)  does  not  repeat  Itself  every  ir/6 
radians  on  the  interval  Cn/6,  511/63,  attention  must  be  directed  to  the 
cycloconverter  connection  scheme  utilized  so  that  the  eliminated  harmonics 
are  not  inadvertently  reintroduced.  Specifically,  a  conduction  path  must 
be  provided  so  that  all  harmonics  remaining  in  the  phase  voltage  can  flow 
if  the  phase  current  is  to  be  of  the  same  harmonic  content.  Two  connection 
realizations  are  possible  wherein  no  restriction  on  phase  current  harmonics 
are  imposed  by  Kirchhoff's  current  law.  The  first  of  these  acceptable 
connections  is  that  of  the  midpoint  cycloconverter  with  the  neutral  lead 
in  place.  A  second  realization  would  be  the  full-bridge  cycloconverter 
with  total  phase  isolation.  Since  this  latter  connection  requires  36  SCRs, 
it  was  not  used  in  the  reported  study  due  to  a  factor  of  three  on  component 
requirement  over  the  midpoint  cycloconverter  case.  However,  it  should  be 
pointed  out  that  the  full-bridge  cycloconverter  with  its  six-pulse  output 
has  a  smaller  ripple  magnitude  and  may  offer  a  slight  efficiency  advantage. 

5.4  Fourier  Spectrum  Analysis 

Nature  of  the  waveforms  for  current  and  instantaneous  developed 
torque  are  such  that  visual  inspection  to  determine  harmonic  content  is 
not  practical  or  accurate.  The  listing  of  a  FORTRAN  program  to  find  the 
Fourier  coefficients  is  given  in  Appendix  D. 

This  program  was  used  to  generate  the  data  for  the  normalized  Fourier 
spectrum  plots  of  Section  III  for  harmonics  up  through  forty-eight  in  order 
to  verify  thap  attempt  to  selectively  eliminate  harmonics  was  successful. 
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SECTION  III 


RESULTS 


1.0  SUMMARY 

Three  speed  conditions  were  selected  for  numerical  study  to  evaluate 
the  effectiveness  of  the  harmonic  elimination  technique  in  removing 
undesired  torque  pulsations  from  the  midpoint  cycloconverter  driven 
brushless  DC  motor  system  when  using  synchronous  envelope  control: 

(a) .  4500  rpm  (10%  speed  case) 

(b) .  2250  rpm  (  5%  speed  case) 

(c) .  450  rpm  (  1%  speed  case) 

Although  these  speeds  are  great  enough  that  elimination  of  harmonics  as 
determined  by  equation  (13)  would  not  be  necessary,  there  are  two  reasons 
for  their  selection.  The  first  is  conservation  of  computer  time.  Computa¬ 
tion  of  a  set  of  data  for  the  450  rpm  (IX  speed  case)  required  approxi¬ 
mately  8  hours  of  CPU  time  on  a  Control  Data  PDF  11/44  computer.  Since 
many  trlal-and-error  runs  were  necessary  for  each  data  set,  runs  at  lower 
motor  speed  become  prohibitive  timewise.  Secondly,  if  the  harmonics  can 
be  successfully  eliminated  at  higher  speeds,  then  there  is  not  difficulty 
to  be  encountered  at  reduced  speed.  Justification  lies  in  the  fact  that 
the  undeslred  harmonics  are  totally  eliminated  from  the  phase  voltage  set. 
Any  re-introduction  of  the  undesired  harmonics  into  the  phase  current 
wave  forms  would  be  attributable  to  commutation  overlap  which  is  dependent 
on  the  La/R  time  constant  that  determines  current  decay  at  the  end  of 
each  phase  current  conduction  sequence.  Since  this  time  constant  is  not 
speed  dependent,  the  commutation  overlap  interval  is  a  larger  portion  of 
the  current  wave  form  period  at  high  speed  than  it  is  at  low  speed.  Thus, 
the  departure  of  phase  current  from  the  form  of  the  phase  voltage  is 
greater  at  high  speed  than  at  low  speed.  It  is  concluded  that  if 
undeslreable  harmonics  are  not  re-introduced  at  high  speed,  then  success 
at  low  speed  is  guaranteed. 

For  each  speed  case,  the  nature  of  current  and  torque  were  examined 
for  application  of  a  balanced  set  of  phase  voltages  of  the  following 
nature: 

(a) .  No  harmonic  eliminated. 

(b) .  Fifth  and  seventh  harmonics  eliminated. 

(c) .  Fifth,  seventh,  eleventh,  and  thirteenth  harmonics  eliminated. 


A  brushless  DC  machine  system  has  an  extra  Input  or  degree  of  freedom 
over  either  a  conventional  DC  machine  or  a  synchronous  machine  drive 
system  in  that  the  angular  displacement  between  the  stator  and  rotor 
magnetic  fields  can  be  specifically  controlled.  For  this  study,  this 
extra  degree  of  control  freedom  is  implemented  by  specification  of  the 
angle  (a)  measured  from  the  onset  of  positive  phase  current  conduction 
to  the  positive  going  zero  crossing  of  the  associated  motor  phase  counter 
emf.  A  convenient  value  (a  -  30°)  is  used  for  most  of  the  study;  however, 
since  selection  of  a  shifts  the  position  of  motor  phase  counter  emf  with 
respect  to  phase  current,  the  nature  of  coomutation  overlap  can  be  altered 
by  the  value  of  a.  The  effect  of  changing  a  over  the  range  from  15°  to 
40°  is  studied  for  one  speed  case. 

2.0  PERFORMANCE  STUDIES 

Equations  (1)  -  (3)  can  be  written  in  compact  matrix  form  as 


-ii 

dt  - 


=-1  CR]  i  +  -1  v*  -  ie 


-  La^ 


La~ 


(20) 


where: 


"*r 

(21)  ;  e  - 

~  P 

sin(u)t-aX 

i  - 

i2 

sin(uit-a-  12L  ) 

3 

_  13 

_sin(cot-a+  ) 

(22) 


-R, 

h(a)t) 

CR]  - 

1 

r2 

(23)  ;  v*  -  vd(t) 

h(u>t  -  ^  ) 

R 

L  3  J 

h(u)t  +  ^  ) 

(24) 


Equation  (20)  is  a  nonlinear  set  of  differential  equations  due  to  the 
elements  of  CR]  being  functions  of  current.  Any  attempt  at  linearization 
for  a  wide  range  of  i_ would  unacceptably  alter  the  physical  nature  of 
the  problem  by  removing  the  rectifying  characteristic  of  the  SCRs.  Conse¬ 
quently,  a  closed  form  solution  of  (20)  is  not  possible  and  a  numerical 
integration  must  be  implemented.  A  fixed  increment  (1.0  x  10~^s)  fourth- 
order  Runge-Kutta  numerical  integration  method  is  utilized  for  this  work. 

After  forming  a  modulation  function  through  use  of  the  programs  of 
Appendix  B,  selecting  a  particular  motor  speed  (Uq  -  2w/p),  and  arbitrarily 
choosing  a  value  of  SCR  firing  angle  ag(see  equation  (17)),  numerical 
solution  to  (20)  is  implemented  subject  to  initial  conditions  i/0)  -  0. 
Integration  is  continued  until  steady-state  is  reached.  Instantaneous 
values  of  developed  torque  are  calculated  according  to  (5)  and  numerically 


averaged  over  the  last  steady-state  cycle  of  integration  to  give  the 
average  value  of  developed  torque.  This  value  of  average  torque  is 
compared  with  a  specified  value  for  the  operating  point.  If  the  average 
value  calculated  is  different  from  the  specified  value,  then  ats  is 
appropriately  adjusted  and  the  numerical  integration  is  repeated  until  the 
difference  magnitude  between  calculated  and  specified  torque  are  suffi¬ 
ciently  close  in  an  epsilon  sense.  For  all  cases  studied,  the  specified 
torque  was  taken  as  1.5  N*m. 

The  listing  of  a  FORTRAN  program  to  compute  the  performance  as 
described  above  is  presented  in  Appendix  D.  In  addition  to  the  discussed  - 
computations,  the  program  also  calculates  the  RMS  value  of  phase  current 

over  the  last  steady-state  cycle  as  (  i.  £  i2  At)*5.  Also,  values  of  phase 

current  and  instantaneous  torque  are  stored  over  the  last  steady-state 
cycle  for  use  by  the  Fourier  Spectrum  Program  of  Appendix  C. 

Results  of  the  various  speed  cases  and  harmonic  elimination  studies 
are  summarized  in  the  paragraphs  that  follow. 

2.1  4500  RPM  (10Z  Speed  Case) 

2.1.1  Unmodulated  Phase  Voltage. 


Input  Data 

Motor  speed,  •  4500  rpm 

Motor  frequency,  f  ■  150  Hz 
SCR  firing  angle,  as  m  58.68° 

Specified  average  torque,  Tjav  -  1.5  N-m 

Output  Data 

Average  torque,  Tdav  ■  1.501  N*m 
RMS  phase  current,  I^g  -  72.77  A 


A  plot  of  phase  current  ij  is  given  by  Figure  6.  The  associated 
Fourier  spectrum  of  11  is  shown  by  Figure  7.  Predominantly  triplen 
harmonic  current  that  flows  through  the  neutral  connection  is  depicted  by 
Figure  8. 

A  plot  of  the  instantaneous  developed  torque  where  the  sixth  harmonic 
is  obvious  is  given  by  Figure  9.  A  Fourier  spectrum  of  the  instantaneous 
torque  is  shown  in  Figure  10  where  existence  of  all  multiples  of  six  as 
theoretically  predicted  are  noted. 


Figure  7.  Fourier  Spectrum  of  Phase  Current  -  4500  RPM,  No  Modulation 


Figure  8.  Neutral  Current  -  4500  RPM,  No  Modulation 
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Figure 


2.1.2  Elimination  of  Sixth  Harmonic  of  Torque.  The  sixth  harmonic  of 
torque  is  eliminated  by  removing  the  fifth  and  seventh  harmonics  of  phase 
current . 


Input  Data 

Motor  speed,  Pq  ■  4500  rpm 

Motor  frequency,  f  -  150  Hz 
SCR  firing  angle,  ct8  -  46.45° 

Specified  average  torque,  Tdav  -  1.5  N-m 
Angles  of  modulation  function,  a.  -  40.76° 

a*  -  47.73° 
o3  -  58.65° 


Output  Data 


Average  torque,  T 
RMS  phase  current 


day 
*  lRMS 


1.5001  N-m 
-  74.24  A 


A  plot  of  phase  current  1^  is  given  by  Figure  11.  The  associated 
Fourier  spectrum  of  1^  is  shown  by  Figure  12,  where  negligibly  small  fifth 
and  seventh  harmonics  are  noted.  Current  flow  through  the  neutral  line 
is  depicted  by  Figure  13. 

A  plot  of  the  instantaneous  developed  torque  is  given  by  Figure  14. 
The  Fourier  spectrum  of  the  developed  torque  is  shown  in  Figure  15  where 
it  is  seen  that  the  sixth  harmonic  is  near  zero. 


2.1.3  Elimination  of  Sixth  and  Twelfth  Harmonics  of  Torque.  Both  the 
sixth  and  twelfth  harmonics  of  torque  are  eliminated  by  removing  the  fifth, 
seventh,  eleventh,  and  thirteenth  harmonics  of  phase  current. 

Input  Data 

Motor  speed,  *  4500  rpm 

Motor  frequency,  f  -  150  Hz 
SCR  firing  angle,  o#  ■  42.63° 

Specified  average  torque,  Tjav  *  1.5  N*m 
Angles  of  modulation  function,  a3  ■  38.73° 

u2  -  42.13° 
a3  -  57.25° 
a4  -  61.93° 
o5  -  66.86° 


Output  Data 

Average  torque,  T<jav  -  1.505  N.m 

RMS  phase  current,  Igjjs  ”  74.59  A 
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Figure  11.  Phase  Current  -  4500  RPM,  Fifth  and  Seventh  Harmonic  Eliminated 


Figure  14.  Developed  Torque  -  4500  RPM,  Sixth  Harmonic  Eliminated 


Sixth  Harmenic  Eliminated 


A  plot  of  phase  current  Is  given  by  Figure  16.  The  associated 
Fourier  spectrum  of  ij  is  shown  by  Figure  17  where  negligibly  small  fifth, 
seventh,  eleventh,  and  thirteenth  harmonics  are  noted.  Current  flow 
through  the  neutral  line  is  depicted  by  Figure  18. 

A  plot  of  the  instantaneous  developed  torque  is  given  by  Figure  19. 
The  Fourier  spectrum  of  the  developed  torque  is  shown  in  Figure  20  where 
it  is  seen  that  the  sixth  and  twelfth  harmonics  are  near  zero. 

2.2  2250  RPM  (5Z  Speed  Case) 


2.2.1  Unmodulated  Phase  Voltage . 


Input  Data 


Motor  speed,  nm  *  2250  rpm 
Motor  frequency,  f  *  75  Hz 
SCR  firing  angle,  8  -  61.93° 

Specified  average  torque,  Tdav  *  1-5  N-m 


Output  Data 


Average  torque,  Tdav  *  1.499  N-m 

RMS  phase  current,  Ipjjs  “  69.53  A 


A  plot  of  phase  current  ij  is  given  by  Figure  21.  The  associated 
Fourier  spectrum  of  i1  is  shown  by  Figure  22. 

A  plot  of  the  instantaneous  developed  torque  is  given  by  Figure  23. 

A  Fourier  spectrum  of  the  developed  torque  is  shown  by  Figure  24. 

2.2.2  Elimination  of  Sixth  Harmonic  of  Torque.  The  sixth  harmonic  of 
torque  is  eliminated  by  removing  the  fifth  and  seventh  harmonics  of  phase 
current. 


Input  Data 


Motor  speed,  nm  ■  2250  rpm 
Motor  frequency,  f  ■  75  Hz 
SCR  firing  angle,  as  -  50.21° 

Specified  average  torque,  Tjav  *  1.5  N-m 
Angles  of  modulation  function,  a.  -  40.76° 

a2  -  47.73° 
a3  -  58.65° 


Output  Data 


Average  torque,  Tda^  *  1.5001  N-m 
RMS  phase  current,  Irms  *  73.96  A 
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Figure  17.  Fourier  Spectrum  of  Phase  Current  -  4500  RPM,  Fifth, 

Seventh,  Eleventh,  and  Thirteenth  Harmonics  Eliminated 


Figure  20.  Fourier  Spectrum  of  Developed  Torque  —  4500  RPMt 
Sixth  and  Twelfth  Harmonics  Eliminated 


Figure  22.  Fourier  Spectrum  of  Phase  Current  -  2250  RPM,  No  Modulation 


250  RPM. 


Fourier  Spectrum  of  Developed  Torque  **  2250  RPM,  Ho  Modulation 


A  plot  of  phase  current  Is  given  by  Figure  25.  The  associated 
Fourier  spectrum  Is  shown  by  Figure  26. 

A  plot  of  the  instantaneous  developed  torque  is  given  by  Figure  27. 
The  Fourier  spectrum  of  the  developed  torque  is  shown  by  Figure  28. 

2.2.3  Elimination  of  Sixth  and  Twelfth  Harmonics  Of  Torque.  Both  the 
sixth  and  twelfth  harmonics  of  torque  are  eliminated  by  removing  the  fifth, 
seventh,  eleventh,  and  thirteenth  harmonics  of  phase  current. 


Input  Data 


Motor  speed,  ■  2250  rpm 

Motor  frequency,  f  ■  75  Hz 
SCR  firing  angle,  aa  *  47.31° 

Specified  average  torque,  T^^  -  1.5  N*m 
Angles  of  modulation  function,  oti  *  38.73° 

a2  -  42.13° 
a3  -  52.25° 
a.  -  61.93° 
a5  -  66.86° 


Output  Data 

Average  torque,  Tdav  ■  1.5003  N*m 
RMS  phase  current,  Igjjs  *  74.64  A 

A  plot  of  phase  current  i^  is  given  by  Figure  29.  The  associated 
Fourier  spectrum  of  1^  is  shown  by  Figure  30. 

A  plot  of  the  instantaneous  developed  torque  is  given  by  Figure  31. 
The  Fourier  spectrum  of  the  developed  torque  is  shown  by  Figure  32. 

2.2.4  Elimination  of  Sixth  and  Twelfth  Harmonics  of  Torque  (a  ■  15°).  Tt 
angle  from  the  onset  of  positive  phase  current  to  the  positive  going 
crossing  of  the  associated  motor  phase  counter  emf  is  changed  from  30°  to 
15°  to  examine  effect  on  current  harmonics. 


Input  Data 


Motor  speed,  n^  ■  2250  rpm 
Motor  frequency,  f  ■  75  Hz 
SCR  firing  angle,  aa  ■  59.84° 

Specified  average  torque,  T<|av  -  1.5  N*m 
Angles  of  modulation  function,  ai  ■  38.73° 

a2  -  42.13° 
a3  -  52.25° 
a4  -  61.93° 
c*5  -  38.73° 
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Figure  25 


Harnonic  Nunber 

Figure  26.  Fourier  Spectrum  of  Phase  Current  -  2250  RPM,  Fifth 
and  Seventh  Harmonics  Eliminated 


Harwonic  NuMber 

Figure  28.  Fourier  Spectrum  of  Developed  Torque  -  2250  RPM 
Sixth  Harmonic  Eliminated 


Degrees 

ie  Current  -  2250  RPM,  Fifth,  Seventh,  Eleventh 
Thirteenth  Harmonics  Eliminated 


Figure  30.  Fourier  Spectrum  of  Phase  Current  -  2250  RPM,  Fifth,  Seventh 
Eleventh,  and  Thirteenth  Harmonics  Eliminated 
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Figure 


Harnonic  Nunber 

Fourier  Spectrum  of  Developed  Torque  -  2250  RPM, 
Sixth  and  Twelfth  Harmonics  Eliminated 


Average  torque,  T<jav  -  1.502  N*m 
RMS  phase  current,  -  53.33  A 


A  plot  of  phase  current  1.  Is  given  by  Figure  33.  The  associated 
Fourier  spectrum  Is  shown  by  Figure  34. 

A  plot  of  Instantaneous  developed  torque  is  given  by  Figure  35. 

The  Fourier  spectrum  of  the  developed  torque  is  shown  by  Figure  36. 

The  re-lntroduced  harmonic  magnitudes  are  reduced  to  a  lower  level 
than  for  a  “  30°. 

2.2.5  Elimination  of  Sixth  and  Twelfth  Harmonics  of  Torque  (a  ■  40°) . 
The  angle  from  onset  of  positive  phase  current  to  the  positive  going 
crossing  of  the  associated  motor  phase  counter  emf  is  set  at  40°  to 
examine  the  effect  on  current  harmonics. 

Input  Data 

Motor  speed,  1%  ■  2250  rpm 

Motor  frequency,  f  -  75  Hz 
SCR  firing  angle,  as  »  16.77° 

Specified  average  torque,  Tjjav  -  1.5  N-m 
Angles  of  modulation  function,  a1  «  38.73° 

a2  -  42.13° 
a3  -  52.25° 

-  61.93° 

a5  -  66.86° 


Output  Data 

Average  torque,  T<jav  -  1.5003  N*m 
RMS  phase  current,  IrmS  "  107.52 


A  plot  of  phase  current  i1  is  given-  by  Figure  37.  The  associated 
Fourier  spectrum  of  phase  current  is  shown  by  Figure  38. 

A  plot  of  instantaneous  developed  torque  is  given  by  Figure  39. 

The  Fourier  spectrum  of  developed  torque  is  shown  by  Figure  40. 

Although  the  angle  power  factor  angle  (a  +  30°)  has  been  extended  to 
75°  resulting  in  a  large  current  to  maintain  the  same  power  flow,  no 
significant  re-introduction  of  undesired  harmonics  is  noted. 
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Harnonic  Huriber 

Fourier  Spectrum  of  Developed  Torque  -  2250  RPM,  a  -  15° 
Sixth  and  Twelfth  Harmonics  Eliminated 


60  120  180  240  390  360 

Degrees 

Figure  37.  Phase  Current  -  2250  RPM,  a  -  40°,  Fifth,  Seventh, 

Eleventh,  and  Thirteenth  Harmonics  Eliminated 
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Fourier  Spectrum  of  Phase  Current  -  2250  RPM,  a  -  40°,  Fifth 
Seventh,  Eleventh,  and  Thirteenth  Harmonics  Eliminated 


Developed  Torque  -  2250  RPM,  a  -  40°,  Sixth 
and  Twelfth  Harmonics  Eliminated 


Harnonic  Nunber 

Fourier  Spectrum  of  Developed  Torque  -  2250  RPM, 
Sixth  and  Twelfth  Harmonics  Eliminated 


2.3  450  RPM  (1%  Speed  Case) 

2.3.1  Unmodulated  Phase  Voltage. 


Motor  speed,  nm  «  450  rpm 

Motor  frequency,  f  *  15  Hz 
SCR  firing  angle,  ag  -  64.52° 

Specified  average  torque,  Tdav  *  1.5  N*m 


Output  Data 


Average  torque,  T^  ■  1.4996  N*m 
RMS  phase  current,  IpMS  “  67.1  A 


A  plot  of  phase  current  ii  is  given  by  Figure  41.  The  associated 
Fourier  spectrum  is  shown  by  Figure  42. 

A  plot  of  Instantaneous  developed  torque  is  given  by  Figure  43.  A 
Fourier  spectrum  of  developed  torque  is  shown  by  Figure  44. 

The  .instantaneous  wave  form  plots  for  the  low  speed  case  do  not  dis¬ 
play  all  of  the  rippled  due  to  the  lack  of  capacity  of  the  plot  routine 
to  store  sufficient  points. 

2.3.2  Elimination  of  Sixth  Harmonic  of  Torque.  The  sixth  harmonic  of 
torque  is  eliminated  by  removing  the  fifth  and  seventh  harmonics  of  phase 
current. 


Input  Data 


Motor  speed,  ng  ■  450  rpm 

Motor  frequency,  f  ■  15  Hz 
SCR  firing  angle,  ag  -  52.53° 
Specified  average  torque,  T^av  ■ 
Angles  of  modulation  function,  a1 

a2 

<*3 


1.5  N-m 

-  40.76° 

-  47.73° 

-  58.65° 


Output  Data 


Average  torque,  T,jav  ■  1.4996  N*m 
RMS  phase  current,  “  75.2  A 


A  plot  of  phase  current  ij  is  given  by  Figure  45.  The  associated 
Fourier  spectrum  is  shown  by  Figure  46. 


Figure  42.  Fourier  Spectrum  of  Phase  Current  -  450  RPM,  No  Modulation 


Har-nonic  Nunber 

Lgure  44.  Fourier  Spectrum  of  Developed  Torque  -  450  RPM,  No  Modulation 
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A  plot  of  instantaneous  developed  torque  is  given  by  Figure  47. 
The  Fourier  spectrum  of  developed  torque  is  shown  by  Figure  48. 


2.3.3  Elimination  of  Sixth  and  Twelfth  Harmonics  of  Torque.  Both  the 
sixth  and  twelfth  harmonics  of  torque  are  eliminated  by  removing  the 
fifth*  seventh*  eleventh*  and  thirteenth  harmonics  of  phase  current. 


Input  Data 


Motor  speed,  nm  *  450  rpm 
Motor  frequency*  f  ■  15  Hz 
SCR  firing  angle,  ag  ■  49.78a 
Specified  average  torque*  Tdav  *  ! 
Angles  of  modulation  function* 

a2 


1.5  N*m 

-  38.73° 

-  42.13° 

-  52.25° 

-  61.93° 

»  66.86° 


Output  Data 


Average  torque*  Td„ 
RMS  phase  current  Ii 


1.507  N>m 
-  77.38  A 


A  plot  of  phase  current  i^  is  given  by  Figure  49.  The  associated 
Fourier  spectrum  of  phase  current  is  shown  by  Figure  50. 


A  plot  of  instantaneous  developed  torque  is  given  by  Figure  51.  The 
Fourier  spectrum  of  developed  torque  is  shown  by  Figure  52. 
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Phase  Current  -  450  RPM,  Fifth  and  Seventh  Harmonics  Eliminated 


Figure  46 


Figure  47 


Harnonic  Nunber 

Fourier  Spectrum  of  Developed  Torque  -  450  RPH 
Sixth  Harmonic  Eliminated 
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SECTION  IV 


DISCUSSION 


1.0  INCREASE  IN  OHMIC  LOSSES 

Since  modulation  of  the  wave  form  definitely  redistributes  the 
harmonic  content  of  the  phase  current  wave  formt  an  examination  of  the 
change  In  winding  ohmic  losses  Is  In  order.  For  the  4500  rpm  case  it 
Is  noted  that  the  RMS  phase  current  is  found  to  be  as  follows: 


(a). 

IRMS  “ 

72.77 

A 

for  no  modulation 

• 

(b). 

IRMS  “ 

74.24 

A 

for  fifth  and  seventh  harmonic  eliminated 

(c). 

XRMS  * 

74.59 

A 

for  fifth,  seventh,  eleventh,  and  thirteenth 
harmonic  eliminated 

The  increase  in  ohmic  losses  assuming  a  constant  winding  resistance  is 
given  by  the  ratio  of  RMS  current  for  a  changed  condition  to  the  RMS 
value  of  current  for  the  unmodulated  wave  form.  It  is  concluded  that 
the  ohmic  losses  are  Increased  by  2 .11  for  the  case  of  elimination  of 
fifth  and  seventh  harmonics.  The  ohmic  losses  are  increased  by  2 .51  for 
the  case  of  elimination  of  the  fifth,  seventh,  eleventh,  and  thirteenth 
harmonics . 

Calculated  data  for  the  2250  rpm  and  450  rpm  cases  show  increases 
of  7  to  15Z;  however,  there  is  no  theoretical  basis  for  greater  Increase 
in  these  cases.  The  reason  for  this  apparent  larger  increase  for  the 
lower  speed  cases  is  a  numerical  problem  due  to  storing  and  calculating 
the  RMS  value  based  on  a  limited  number  of  data  points  (2000)  for  trans¬ 
mittal  to  the  plot  routine.  As  a  result,  the  high  frequency  ripple 
currents  are  sampled  at  irregular  intervals  giving  rise  to  Inaccuracy  in 
calculation  of  the  RMS  value.  For  the  case  of  4500  rpm.  the  data  point 
limit  allows  storage  of  several  points  during  each  ripple  cycle,  thus,  the 
resulting  RMS  values  are  considerably  more  accurate  than  for  the  other 
two  speed  cases. 

2.0  CYCLOCONVERTER  MODE  CHANGE 

As  discussed  and  illustrated  by  Figure  2,  a  neutral  connection  is 
necessary  to  successfully  eliminate  selected  phase  current  harmonics. 

Since  harmonic  elimination  is  only  necessary  at  speeds  near  zero,  and 
since  an  ohmic  loss  penalty  is  incurred  as  a  result,  it  might  be 
deslreable  to  not  operate  with  PWM  control  over  the  full  range  of  speed. 
The  Harmonic  Eliminator  Switch  in  the  neutral  line  of  Figure  2  could  be 
opened  for  higher  speeds  and  the  cycloconverter  mode  of  operation  changed 
to  phase-control  from  synchronous  envelope. 
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3.0  EXTENSION  TO  HIGHER  HARMONIC  ELIMINATION 

For  specific  illustration,  this  study  has  only  considered  elimination 
of  torque  harmonics  through  the  twelfth.  As  motor  speed  approaches  zero, 
it  may  be  necessary  to  eliminate  even  higher  torque  harmonics.  Theoreti¬ 
cally,  there  is  no  limitation  to  be  encountered. 


SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  work  performed  by  this  study  shows  that  systematic  reduction  of 
harmonic  torque  pulsations  in  brushless  DC  drives  is  feasible  by  use  of 
selective  current  harmonic  elimination.  Further,  there  appears  to  be 
only  a  small  ohmic  loss  penalty. 

Further  theoretical  study  should  be  made  to  properly  assess  quanti¬ 
tatively  the  effect  of  increased  motor  leakage  inductance  and  of  interphase 
reactors  on  "reintroducing"  the  harmonics  eliminated  from  the  phase 
voltage  due  to  extending  the  commutation  overlap. 

Also,  a  hardware  realization  of  a  midpoint  cycloconverter  brushless 
DC  drive  using  the  control  principles  formulated  in  this  study  should  be 
made  to  verify  the  results. 
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APPENDIX  A 


HARMONICS  OF  PWM  TIME  FUNCTION 


A  modulation  function  h(oit)  was  derived  in  Appendix  B  and  used  in 
Section  II  -  5.0  to  form  a  modulated  phase  voltage.  For  example. 

Van  ■  h(ut)  vd(t). 

Harmonic  content  of  a  modulated  voltage  is  of  concern.  The  study 
is  best  understood  by  use  of  a  typical  example.  Let  h(wt)  contain  all 
odd  harmonics  except  the  fifth  and  seventh,  then 


h(uit)  *  a^cos  ait  +  a^cos  3wt  +  agcos  9wt  +  a^cos  Hut  +  ••• 

+  ancos  nut  +  •••  (A.l) 

Assume  the  vd(t)  contains  a  DC  term  plus  one  high  frequency  harmonic,  then 


vd(t)  ■  bQ  +  bmcos  mut  (A. 2) 

Forming  v^  as  the  product  of  h(ut)  and  v^Ct)  and  simplifying  yields 

vd(t)  ■  Cjcos  ut  +  c3co8  3ut  +  c9cos  9ut  +  c^cos  llut  +  ••• 

+  d^jcos  (m+l)ojt  +  dTn_^cos  (m-l)ut  +  d^)  ^cos  (nH-3)ut 
+  dm_3COS  (m-3)wt  +  •  •  •  (A.  3) 


It  is  seen  from  (A. 3)  that  through  the  cross  product  terms,  the  PUM 
expression  may  contain  frequencies  that  were  selectively  eliminated  from 
h(o)t).  However,  the  coefficients  (d^fi,  dm_i,  ...)  of  (A. 3)  are  products 
of  Sq  and  bg,.  Since  decreases  as  l/n,  if  m  is  large,  then  the 
coefficients  of  the  low  frequency  components  of  h(ut)  resulting  cross 
product  terms  will  be  negligibly  small.  It  is  concluded  that  the  fre¬ 
quency  of  the  AC  source  must  be  large  compared  with  the  motor  frequency 
if  the  PWM  function  is  to  not  contain  the  harmonics  selectively  eliminated 
from  the  modulating  function  h(ut).  But,  at  near  zero  speed,  the  source 
frequency  is  large  compared  with  the  motor  frequency. 
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APPENDIX  B 


MODULATION  FUNCTION  PROGRAMS 


1.0  INITIAL  SOLUTION 

As  first  procedure,  a  candidate  or  Initial  modulation  function,  h(oat) , 
Is  found  to  serve  as  a  basis  for  an  optimization  search  to  determine  the 
final  modulation  function. 

1.1  Theory  of  Harmonic  Elimination 

Extending  Patel  and  Hoft's  pulse  width  modulation  technique  C14D  to 
Include  a  slack  variable,  it  is  possible  to  find  a  modulation  function 
that  has  selected  harmonics  elimination  while  existing  on  the  quarter  wave 
interval  from  30°  to  90°. 

The  Fourier  series  representation  of  the  modulation  function  assuming 
odd  quarter  wave  symmetry  is 


f(ut)  «  £  ansin(nut)  ,  for  odd  n 

1 

a2  a4  aM 

/„  sin(nut)dut  +  f  sin(nwt)du>t  +  •••  L  sin  (nut)  dm  t,  (B.l) 
1  “3  Ti-1 


30°  <  ax  <  a2  <  ...  <  <  90° 


This  reduces  to 


_4 


I  (-l)k+1cos(nak) 
k-1 


(B.2) 


for  M  both  odd  and  even.  The  number  of  harmonics  eliminated  is  in  general 
the  same  as  the  number  of  switching  transitions  per  quarter  cycle.  For  M 
greater  than  2  the  equations  cannot  be  solved  in  closed  form,  so  numerical 
methods  are  needed  to  obtain  solutions.  The  following  normalized  system 
of  equations  is  to  be  used  to  obtain  solutions: 


fnl  "  I  (-l)k+1cos(nlCik)  -  o 

M  ... 

fn2  "  I  (-D  cos(n2ak)  -  0 

:  m 

fnM  "  l  (-D^cosCtyijj)  •  0 


(B.3) 


where  n^,^, •**nM  are  the  harmonics  to  be  eliminated. 

Many  of  the  solutions  of  (B.3)  are  meaningless  or  impractical,  so  it 
is  necessary  to  search  by  brute  force  using  a  loose  convergence  criterion 
for  candidate  solutions,  and  then  to  use  a  numerical  method  such  as  Newton 
Raphson  to  obtain  a  refined  solution. 

The  basic  method  of  selecting  candidate  modulation  function  for 
further  optimization  is  illustrated  by  the  flow  chart  of  Figure  B.l.  A 
candidate  set  is  one  that  meets  the  following  criteria: 


^squar  ^  fmin 

f squar  ”  fnl  +  fn2  ‘ *  *  fnM 

If  for  the  candidate  modulation  function,  30°  <  aj  <  a2  <  *  *  *  <  aM  <  90°, 
then  the  solution  is  considered  for  optimization.  Otherwise,  the  solution 
is  discarded. 


Figure  B.l.  Flow  Chart  for  Selection  of  Candidate 
Modulation  Functions 


>iV*v 


.  /L 


1.2  Initial  Solution  Program 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  main  program  c 

c  c 

ccccceeccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
on-  number  of  harmonics  to  be  eliminated 
c  this  includes  the  slack  variable 

c  inc  -  this  is  the  amount  that  each  alpha 

c  will  be  Incremented  by  during  the 

c  infinite  search 

c  fmin  -  newton  rapson  method  is  used  to 

c  search  for  a  solution  if  the  sum 

c  of  squares  is  less  than  fmin 

c  r(l)  -  these  are  the  harmonics  to  be  reduced 

c  amin  -  this  is  the  minimum  value  of  alpha 

c  for  the  infinite  search 

c  iwr  -  not  used 

c  del  -  not  used 

dimension  a(10)yr(10),f(10) 

common  plfconv,iwr,del 

real  inc 

call  sopen(O.O) 

pi*3. 141592653 

klminaO 

k2min=0 

k3min=0 

kUminsO 

k5mins0 

k6min=0 

kl =1 

k2*1 

k3*1 

k4sl 

k5=1 

k6*1 

conv* 180. 0/pi 

read(9,*)  n 

write(6#*)  1  ns’,n 

read(9,*)  inc 

write (6,*)  *  inea'.inc 

read(9,*)  fmin 

write(6,*)  1  fmina',fmin 

read(9»*)  (r(i),isi,7) 

write(6,#)  *  harmonicsa*t(r(i),ial,n) 

read (9 »•)  amin 

write(6t*)  '  amina»,amin 

read(9,*)  iwr 

write(6,»)  *  iwra’,iwr 


max=int(  (90.0-amin)/inc  ) 

k1max=lnt(  (90.0-amln)/lnc  ) 

k2naxsmax 

k3max=max 

kUnaxsmax 

k5maxsmax 

k6naxsoax 

isep=lnt(  1.0/lnc  )  +1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

e 

c  the  following  are  essentially  do  loops  but  were 

c  necessary  to  get  past  a  problem  with  the  compiler 

c 

klminsl 

100  k1=k1min-1 

105  k1=k1+1 

lf(  n  .eq.  1  )  go  to  1000 

k2mln&k1+lsep 
200  k2=k2min-1 

205  k2*k2+1 

lf(  n  .eq.  2  )  go  to  1000 

k3mln=k2+isep 
300  k3-k3mln-1 

305  k3*k3+1 

lf(  n  .eq.  3  )  go  to  1000 

k4mln=k3+lsep 
400  k4=k4min-1 

405  k4sk4+1 

lf(  n  .eq.  4  )  go  to  1000 

k5minak4+laep 
500  k5«k5mln>1 

505  k5*k5+1 

lf(  n  .eq.  5  )  go  to  1000 

k6mln«k5+lsep 
600  k6sk6mln-1 

606  k6sk6+1 

If C  n  .eq.  6  )  go  to  1000 

1000  a(1)*(  float (k 1 )*inc+amin  )/conv 

a(2)«(  float (k2)#inc+amln  )/oonv 
a(3)*(  float(k3)*inc+amin  )/conv 
a(4)»(  float (k4)*lnc-«>amln  )/conv 
a(5)«(  float(k5)*inc+amin  }/conv 
a(6)«(  float (k6)*lnc-»-amln  )/conv 
call  ff(  a,r,f,n  ) 
oall  square (  ftfsquar,n  ) 

If C  fsquar  .It.  fmln  )  call  look(  a,r,f ,fsquar,n  ) 
lf(  k6  .It.  k6max  .and.  n  .ge.  6  )  go  to  606 

lf(  k5  .It.  k5max  .and.  n  .ge.  5  >  go  to  505 

lf(  k4  .It.  k4max  .and.  n  .ge.  4  )  go  to  405 

lf(  k3  .It.  k3max  .and.  n  .ge.  3  )  go  to  305 


if(  k2  .It.  k2max  .and.  n  .ge.  2  )  go  to  205 
if(  kl  .It.  klmax  .and.  n  .ge.  1  )  go  to  105 
wrlte(6, 3100) 

write (6, 3000)  (  a(i)*conv,i=1,n  ) 
call  sclose(O.O) 

3000  format (  10(e15.7,5x)  ) 

3100  formatC *  a*') 

3200  formatC  harmonicas') 

stop 
end 
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ccoccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccc 
cocccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  a  c 

c  c 

ccccccccceccececccccoccocccccccccccceccccccccccccccccccccccccccccccecccc 
cccccccccoceococccceccccccoccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  fa(  a,da,n  ) 
dimension  a(10),da(10) 
common  pi»conv,iwr,del 

do  200  ls1,n 
a(i)=a(i)-da(i) 

200  continue 
return 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  close  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  sclose(dummy) 
closed ) 
close(2) 
close (3) 
closed) 
close (7) 
close(8) 
closed) 
return 
end 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  C 

c  subroutine  df  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  fdf(  a,r,f,df,n  ) 
dimension  a(10),r(10),f(10),df(10,10) 
common  pi,conv,iwr,del 


do  200  i*1,n 
signs  1.0 

do  100  Jsl.n 

df(i, j)s-sign*r(i)*sin(  r(i)*a(j)  ) 
signs-sign 

100  continue 

200  continue 
return 
end 


cccccccccecccccccccccccccccccccccccccccccccccocccccccccccccccccccccccccc 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  e 

c  subroutine  f  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccocccc 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccc 
subroutine  ff(  a,r,f,n  ) 
dimension  a(10),r(10),f(10) 
common  pi,convtiwr,del 
do  200  i=1,n 
signsl.O 
f (i)=0.0 

do  100  j=1tn 

f(i)sf(i)+sign*cos(  r(i)»a(j)  ) 
slgns-sign 

100  continue 

200  continue 


return 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  look  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  look(  a,r,f .fsquar ,n  ) 

dimension  da(10),r(10),a(10),f(10),df(10,10),dfi(10,10) 

common  pi,conv,iwr,del 

write (6, 3050) 

write(6,3050) 

write(6,3100) 

write(6,3000)  (  a{i)*conv,i=1,n  ) 
write(6,3150) 

write (6, 3000)  (  f(i),i=1,n  ) 

write(6,3200) 

write (6 ,3000)  fsquar 

kount=30 

fminaO.01 

delta=0.001 


do  2500  loops l.kount 
call  fdf(  a,r,f,dftn  ) 
call  matv(  df,n,det,dfi  ) 
call  rnultK  f,dfi,da,ntn  ) 
call  fa(  a,da,n  ) 
call  ff(  a,r,f,n  ) 
call  square (  f, fsquar ,n  ) 

if(  loop  .gt.  5  .and.  fsquar  .It.  fmin  )  go  to  2900 
2500  continue 

2900  write (6, 3050) 

write(6,3100) 

write(6,3000)  (  a(i)*conv,is1,n  ) 
write(6,3150) 

write (6 ,3000)  (  f(i),is1,n  ) 
write (6, 3200) 
write (6 ,30000  fsquar 
3000  format(  10(e10.4,5x)  ) 

3050  format (  ’  ’  ) 

3100  format(  ’  a= '  ) 

3150  format(  *  fs1  ) 

3200  format(  1  fiquars*  ) 

4000  format ( ’test’ ,13) 

5000  return 

end 


2/2 


AD-A144  216 
UNCLASSIFIED 


HARMONIC  CONTROL  TO  REDUCE  TORQUE  PULSATIONS  IN 
BRUSHLESS  DC  MOTOR  DRIVES(U)  KENTUCKY  UNIV  RESEARCH 
FOUNDATION  LEXINGTON  J  J  CATHEY  MAR  84 
AFOSR-TR-84-060?  AFOSR-82-8189 


NL 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  SUBROUTINE  MATV  C 

c  c 

c  c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MATV(A,N,DET,B) 

DIMENSION  INDEX (20,2), IPVOT <20>,A(10,10),B(10,10), PIVOT ( 20 ) 

ONES  1.0 
ZEROsO.O 


c 

DO  10  1=1, N 

DO  5  J=1,N 
B(I,J)=A(I,J) 

5  CONTINUE 

10  CONTINUE 
DETsONE 
DO  15  J=1,N 
15  IPVOT( J)=0 

DO  75  1=1, N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  FOLLOWING  15  STATEMENTS  C 
C  FOR  SEARCH  FOR  PIVOT  ELEMENTS  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TsZERO 
DO  40  J=1 , N 

IF(IPV0T(J)-1)  20,40,20 
20  DO  35  Ksl.N 

IF(IPVOTCK)-I)  25,35,100 
25  CONTINUE 

DUABS(T) 

D2=ABS(  B(J,K)  ) 

IF(D1-D2)30,35,35 
30  IROWsJ 

ICOLsK 
T*B(J,K) 

35  CONTINUE 

40  CONTINUE 
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IPVOT (ICOL )=IPVOT (ICOL )+1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  FOLLOWING  10  STATEMENTS  PUT  C 
C  PIVOT  ELEMENT  ON  DIAGONAL  C 
C  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF(IROW-ICOL)  45,55,45 
45  DETs-DET 

DO  50  L=1 , N 
TsB(IROW,L) 

B(IROW,L)sB(ICOL,L) 

50  B(ICOL,L)=T 

55  INDEX(I , 1 )=IROW 

INDEX(I,2)=IC0L 
PIVOT(I )=B(ICOL , ICOL ) 

DET  sDET*PIVOT (I ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FOLLOWING  3  STATEMENTS  TO  DIVIDE  PIVOT  C 

C  ROW  BY  PIVOT  ELEMENT  C 

C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

B ( ICOL , ICOL ) =ONE 
DO  60  L=1,N 

60  B (ICOL, L )sB (ICOL, L) /PIVOT (I ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"vCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FOLLOWING  7  STATEMENTS  TO  REDUCE  C 
C  NON-PIVOT  ROWS  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DO  75  LIsl.N 
IF (LI -ICOL)  65,75,65 
65  T*B(LI,ICOL) 

B ( LI, ICOL )= ZERO 
DO  70  L=1,N 

70  B(LI,L)sB(LI,L)-B(ICOL,L)»T 

75  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FOLLOWING  11  STATEMENTS  TO  C 
C  INTERCHANGE  COLUMNS  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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g.y  ,»  ITOffgTOWBWWW  mJWHBJHfffiyvgoyW’iV  •v'.vj*v*.'.n'w.wmwx 


DO  95  1=1, N 
LsH-I+1 

IF ( INDEX (L, 1 )-INDEX(L,2)  )  85,95,85 
85  JROWsINDEX(L.I) 

JC0L=INDEX(L,2) 

DO  90  Ks1,N 

T=B(K, JROW) 

B(K, JR0W)=B(K, JCOL) 

B(K, JC0L)=T 
90  CONTINUE 

95  CONTINUE 

100  RETURN 

END 


'.i  1 1 *r  r ,.-*.vi.t->i-t'.  ■»;,■*.■» . '  :■/ 
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cccceccccccccccccccccccccecccccccecccccccccccccccccccccccccccccccccccccc 
cceccccccccecccccccccccccccccccccccccccccccccccccccccccccccccceccccccccc 
c  c 

c  subroutine  multiply  c 

c  c 

cccccccccecccccccccececccocccccceccccccccccecccccccccccceccccccccccccccc 
cceccccccccccccccccccecccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  multi (  a,b,c,10,k0  ) 
dimension  a(10),b(10, 10),e(10) 
do  200  is 1,10 
c(i)s0.0 
do  100  ks1,k0 

c(i)sc(i)+b(i,k)*a(k) 


100  continue 

200  continue 
return 
end 
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cceccccecccccceccecceccccecccccccecccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  open  c 

e  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccceecccceccc 
cccccccecccccecccccccccccccccccccccccccccccccccccecccccccccecccccccccccc 
subroutine  sopen( dummy) 
open ( unit =9, files "input") 
rewind (9) 
return 
end 


cecccecceccccccecccceecccccccccccccccccccccccccccccccccccccccccccccccccc 
eccceoecececcococeccoecccccccccccceccccccccccccccccccccccccccccccccccccc 
e  c 

c  subroutine  square  c 

o  c 

cccocceeeccccccceceecccccccccceececccccocccccccccceccccccceccccccecccccc 
ccccccceeeececeeoooooeeccoeecocecccccccoeccecccecccccccccccccccccccccccc 
subroutine  squareC  f.fsquar.n  ) 
dimension  f(10) 
fsquarsO.O 
do  200  isl.n 

fsquarafsquar+f(i)*f(i) 

200  continue 
return 
end 


2.0  OPTIMIZATION  OF  MODULATION  FUNCTION 


2.1  Discussion  of  Procedure 

A  candidate  solution  for  the  modulation  function  is  selected  from 
the  output  of  the  Initial  Solution  Program  for  optimization.  The  slack 
harmonic  of  the  Initial  solution  is  changed  in  small  increments  and  fed 
to  a  Newton-Raphson  routine  to  find  the  new  set  of  solutions.  A  con¬ 
tinuous  set  of  o's  results  from  this  procedure  having  "nice"  properties. 
First,  the  new  set  of  solutions  is  well  behaved  in  that  meaningless  sets 
of  solutions  are  not  encountered  and  convergence  is  quick.  Second,  the 
procedure  gives  a  "visual"  picture  of  the  behavior  of  the  solutions  in 
relation  to  the  slack  variable. 

Since  infinite  number  of  modulation  functions  exist,  some  method  is 
needed  to  select  the  "best"  solution.  The  criteria  used  in  this  study 
was  minimizing  the  sum  of  the  squares  of  the  magnitudes  of  the  next  set 
of  current  harmonics  to  be  eliminated  beyond  those  of  present  concern. 
For  example,  if  the  fifth  and  seventh  harmonics  were  being  eliminated 
from  h(  t),  then  the  procedure  is  to  minimize  the  sum  of  the  square  of 
the  magnitudes  of  the  eleventh  and  thirteenth  harmonics. 

A  flow  chart  illustrating  optimization  of  the  modulation  function 
is  given  by  Figure  B.2. 

Figure  B.3  illustrates  the  fora  of  the  modulation  function  for  the 
case  of  elimination  of  fifth  and  seventh  harmonics  and  for  the  case  of 
elimination  of  fifth,  seventh,  eleventh,  and  thirteenth  harmonics. 
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CANDIDATE  SOLUTION 


Figure  B.2.  Flow  Chart  for  Optimizing 
Modulation  Function 
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Figure  B.3.  Form  of  Modulation  Functions. 

(a)  Fifth  and  Seventh  Harmonics  Eliminated 

(b)  Fifth,  Seventh,  Eleventh,  and  Thirteenth 

Harmonics  Eliminated 
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2.2  Modulation  Function  Optimization  Program 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
eccceccecceccccccccccccccccccccceccecccccccccccccccccccceccccccccccccccc 
c  c 

c  main  program  c 

c  c 

cccccccccccccccccccccocccccceccccccccccceccccccccccecccccccccccccccccccc 
ccoccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccc 
c  n  number  of  harmonics  to  be  eliminated 

c  n+1  harmonics  are  eliminated 

c  r(i)  -  harmonics  to  be  eliminated 

c  a(i)  -  these  are  the  angles  from  the  infinite 

c  search  .i.e.  this  is  the  seed 

c  rmax  -  upper  limit  of  slack  variable 

c  rmin  -  lower  limit  of  slack  variable 

c  points  -  number  of  solutions  or  points  to  be  saved 

c  msq  -  number  of  harmonics  to  be  squared  and  sumed 

c  tr(i)  -  harmonics  to  be  summed  and  squared  excluding  tr(1) 

c  tr(1)  -  any  misc.  harmonic  that  the  user  wishes  to  view 

c  incre  -  used  for  increasing  or  decreasing  the  slack  variable 

c  incresl  Increment  the  slack  variable  from  min  to  max 

c  IncresO  decrement  the  slfcck  variable  from  max  to  min 

c  hmag  -  this  is  the  sum  of  sqaures 

dimension  a(10)#r(10)#f(10),h(10),tr(10) 

common  pltconv 

call  sopen(O.O) 

pi* 3. 1^1592653 

conv* 180. 0/pi 

read(9,*)  n 

write(6,»)  '  ns*fn 

read(9,*)  (r(i)fislf5) 

write(6,*)  '  harmonicas', (r(i)» is 1,n) 

read(9t*)  (a(i),isl,5) 

write(6,B)  '  alphass' ,(a(i),is1,n) 

read(9»*)  rmax 

write(6,*)  '  rmax =', rmax 

read (9,*)  rmin 

write(6,*)  *  rmins'trmin 

read(9*B)  points  * 

wrlte(6,*)  '  pointss' .points 

read(9,#)  msq 

write(6,*)  '  msqs',msq 

read(9,»)  (tr(i),isif5) 

write(6,*)  '  hsquar  harmonicas' , (tr(i) ,is2,msq+1 ) 
read(9t*)  incre 

write(6,B)  '  lncres' , incre 


rinos(  rmax-rmin  )/points 
maxrsint(  rmax/rinc  ) 
minrsintC  rmin/rinc  ) 
irmaxsmaxr-minr 
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do  50  i=1 ,5 
a(i)sa(i)/conv 
50  continue 

do  100  irsminr,maxr 

if(  incre  .eq.  1  )  r(n+1 )= float (ir)»rine 
if(  Incre  .eq.  0  )  r(n+1 )=float(  maxr+minr-ir  )*rinc 
call  ff(  a,r,f,n+1  ) 
call  square (  f , fsquar,n+1  ) 
call  look(  a,r,f ,fsquar,n+1  ) 
call  shar(  a,tr,h,n+1 ,msq+1  ) 
call  shsq(  h,hsquar,msq  ) 
hoag=sqrt(  hsquar  ) 
writeO,*)  (  a(i)#conv,isl tn+1 ) 
write(2, 1000)  h(1) 
vrrite(3, 1000)  h(2) 
write(4, 1000)  h(3) 
write(7,1000)  r(n+1) 
write(8, 1000)  hoag 
100  continue 

call  sclose(O.O) 

1000  format(  e10.4  ) 

3000  format (  10(e15.7,5x)  ) 

3100  formate  a*') 

3200  formate  harmonicas') 

stop 
end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  a  c 

c  c 

ccccccccccccccccccccccccccccccccccccccoccccccccccccccccccccccccccccccccc 
cccccccccccecccccceccoccccccccecccccecoccccccccccccccccccccccccccccccccc 
subroutine  fa(  a,da,n  ) 
dimension  a(10),da(10) 
common  pi.conv 

do  200  1=1, n 
a(i)=a(i)-da(i) 

200  continue 
return 
end 
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cccccccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccc 
cccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  close  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  sclose(dunray) 
closed) 
close(2) 
close (3) 

closed) 
close (7) 
closed) 
close (9) 
return 
end 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  C 

c  subroutine  df  c 

C  C 

cccccccccccocccccccccecccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  fdf(  a,r,f,df,n  ) 
dimension  a(10),r(10),f(10),df(10, 10) 
common  pi,conv 


do  200  is1(n 
signs  1.0 

do  100  J=1,n 

df(i, j)s-sign*r(i)*sin(  r(i)»a(J)  ) 
slgnx-sign 

100  continue 

200  continue 
return 
end 
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cccccceccccececccccccccccccccccccccccccccccccccccccccccccccccccccceccccc 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


subroutine  f 


c  c 

ecceeccccoccccccccccccccccccccccccccccoccccccccccccccccccccccccccccccccc 


cccccccceccccccccccccccccccccccccecccccccccccccccccccccccccccccccccccccc 
subroutine  ff(  a,r,f,n  ) 
dimension  a(10),r(10),f(10) 
common  pitconv 
do  200  islyn 
signal. 0 
f ( i)=0. 0 

do  100  jal.n 

f(i)sf(i)+sign*cos(  r(i)*a(j)  ) 
signx-sign 

100  continue 

200  continue 
return 
end 


iv-* 


NSW 


TO 


%s\ 


I 

Sl\| 


m 


.y, 

\ 

,  >i 


it 


SI 


SI 


s 


1 


$ 


cccccccccccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccc 
cceccccccccccccccccccccccccccccccccceecccccccccccecccccccccccccccccccccc 
c  c 

c  subroutine  harmonics  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  shar(  a,tr,h,n,m  ) 
dimension  a(10),h(10),tr(10) 
common  pi,conv 


do  200  is1tm 
signs  1.0 
h(i)s0.0 

do  100  jsl.n 

hdJshdj-tJI.O/f  tr(i)*pi  )*sign»cos(  trd)*a(j)  ) 
slgns-sign 
continue 

continue 

return 

end 
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cccccccccccccccccecccccceccccccececeececccccccccccccccccccccccccccoccccc 
ccccccccccccccccccccccccccccceccccecccoccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  hsquare  c 

c  c 

ccccccccccccccceccccececcccccececccccccecccccccccccceccccccccccecccccccc 
ccccccccccccccccccccccccccceeocceecccccceccccccccccceccccccccccccccccccc 
subroutine  shsq(  h,hsquar,n  ) 
dimension  h(10) 
hsquarsO.O 
do  200  is2tn+1 

hsquar=hsquar+h(i)*h(i) 

200  continue 
return 
end 


ccececcccccccecccecccccccccccccccccccecccecccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
o  c 

c  subroutine  look  c 

e  c 

ecccoecccceccccccccccccccccccccceecceccccccccccccccccccccecccccceccccccc 
ccecoccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccc 
subroutine  look(  a,r ,f , fsquar ,n  ) 
dimension  a(10),r(10),f(10)#df(10, 10),dfi(10, 10) 
common  pi, con v 
kount=30 
fmin*0.01 

do  2500  loops l.kount 
call  fdf(  a,r,f,df,n  ) 
call  matv(  df,n,det,dfl  ) 
oall  mult1(  f,dfl,da,n(n  ) 
call  fa(  a,da,n  ) 
call  ff(  a,r,f,n  ) 
call  square (  f,fsquar,n  ) 

If (  loop  .gt.  5  .and.  fsquar  .It.  fmln  )  go  to  2900 
2500  continue 

2900  dummy =0.0 

3000  format (  10(e15.7,5x)  ) 

3050  formate  'O'  ) 

3100  formate  '  as'  ) 

3150  formate  '  fs'  ) 

3200  formate  '  fsquar* ’  ) 

4000  formate 'test' ,13) 

return 
end 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  MATV  C 

c  c 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MATV(A,N,DET,B) 

DIMENSION  INDEX (20,2), IPVOT ( 20 ) , A ( 1 0 , 1 0 ) , B ( 1 0 , 1 0 ) , PI VOT ( 20 ) 

ONES 1.0 
ZEROsO.O 


DO  10  1=1, N 


DO  5  J= 1 , N 
B(I,J)=A(I,J) 

5  CONTINUE 

10  CONTINUE 

DETsONE 
DO  15  Jsl.N 
15  IPVOT(J)=0 

DO  75  1=1, N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  FOLLOWING  15  STATEMENTS  C 

C  FOR  SEARCH  FOR  PIVOT  ELEMENTS  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

TsZERO 
DO  40  J=1 , N 

IF ( IPVOT (J)-1)  20,40,20 
20  DO  35  Ks1,N 

IF(IPV0T(K)-1)  25,35,100 
25  CONTINUE 

DlsABS(T) 

D2=ABS(  B( J,K)  ) 

IF(D1-D2>30,35,35 
30  IROWsJ 

ICOLsK 
T=B(J,K) 

CONTINUE 


35 


40 


CONTINUE 
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IPVOT ( ICOL ) * IPVOT C ICOL ) +1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  FOLLOWING  10  STATEMENTS  PUT  C 
C  PIVOT  ELEMENT  ON  DIAGONAL  C 
C  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
IFCIROW-ICOL)  45,55,45 
45  DET*-DET 

DO  50  La1,N 
T>B(IROW,L) 

B(IROW,L)aB(ICOL,L) 

50  BCICOL,L)sT 

55  INDEX(I, 1 )sIROW 

INDEX(I,2)«IC0L 
PIVOT ( I )*B( ICOL , ICOL ) 

DETaDET*PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  FOLLOWING  3  STATEMENTS  TO  DIVIDE  PIVOT  C 

C  ROW  BY  PIVOT  ELEMENT  C 

C  C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

B(ICOL, ICOL )= ONE 
DO  60  La1,N 

60  B(ICOL,L)aB(  ICOL,  D/PIVOT  (I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  FOLLOWING  7  STATEMENTS  TO  REDUCE  C 
C  NON-PIVOT  ROWS  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DO  75  Llal.N 
IF (LI -ICOL)  65,75,65 
65  TaB(LI.ICOL) 

B (LI, ICOL) a ZERO 
DO  70  La1(N 

70  B(LI,L)aB(LI,L)-B(ICOL,L)*T 

75  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FOLLOWING  11  STATEMENTS  TO  C 
C  INTERCHANGE  COLUMNS  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


DO  95  1*1, N 
LsH-I+1 

IF(INDEX(L, 1 )-INDEX(L,2)  )  85,95,85 
85  JROW*INDEX(L, 1 ) 

JCOL=INDEX(L, 2) 

DO  90  K=1,N 

TsB(K.JROtf) 

B(K, JROW)=B(K, JCOL) 

B(K, JCOL)=T 
90  CONTINUE 

95  CONTINUE 

100  RETURN 

END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccceccccccccccccccccccccccccccccccecccccccceccccccccccccccccccccc 
c  c 

c  subroutine  multiply  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccceccccccecccc 
eccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 
subroutine  multi (  afbtc,10,k0  ) 
dimension  a(10),b(10, 10),e<10) 
do  200  isl.iO 
c(i)s0.0 
do  100  k=1,k0 

c(i)sc(i)+b(i,k)*a(k) 

100  continue 

200  continue 
return 
end 


■  ■< 
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eccccoccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccecccccoccccccccccccccccccccccccoccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  open  c 

c  c 

ccccccccccocccoccccoccccecccoccccccecccccccccccccccccccccccccccccccccccc 
cccececccccceeccccoocccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  sopen( dummy) 
open(unit=1 ,file="dataO/alpha") 
rewind (1 ) 

open(units2, f ile="dataO/fundamental" ) 
rewind (2) 

open(units3,files"data0/har2tt) 
rewind (3) 

open (units4, f ile="dataO/har 3" ) 
rewind (4) 

open (units7 , f i le="dataO/imghar " ) 
rewind (7) 

open (unit=8 >  f ile= "dataO/hmag" ) 
rewind (8) 

open ( un i t =9 » f i les " input tt ) 

rewind (9) 

return 

end 
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cccccccccccccccccccccccccccccccocccecccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccecccccccccccccccceecccccccccccccccccc 
c  c 

c  subroutine  square  c 

c  c 

ccccccccccccceecccccccceccccccceccccccececcccccccccccccccccceccccccccccc 
cccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccccceccccccc 
subroutine  square (  f,fsquar,n  ) 
dimension  f(10) 
fsquarsO.O 
do  200  i=1,n 

fsquarsfsquar+f(i)*f(i) 

200  continue 
return 
end 
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occccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
0  c 

c  subroutine  square  c 

c  c 

cccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  square(  f,fsquar,n  ) 
dimension  f(10) 
fsquarsO.O 
do  200  i=1,n 

f squar s f aquar +f ( i ) * f ( i ) 

200  continue 
return 
end 
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APPENDIX  C 


FOURIER  SPECTRUM  PROGRAM 


cccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

c 

c 


tor 

X 


mag 

n 

inv 


frequency  spectrum 


c 

c 

c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

-  input  data  (this  is  a  real  variable) 

-  complex  variable  containing  the  input  data 
before  the  fft  subroutine  call  and 
the  complex  frequency  after  the  fft  subroutine  call 

-  magnitude  of  the  complex  frequency 

-  number  of  data  points  (must  be  a  power  of  2) 

-  used  by  fft 

inv=0  fft 

inv=1  inverse  fft 

maxi  -  used  to  decided  highest  frequency  to  be  plotted 

avemag  -  used  to  normalize  data  to  100%  of  the  largest 
harmonic 


dimension  tor(1030) 
complex  x(1030),cmplx 
real  mag 
call  sopen(O.O) 


pi*3. 1415926535 

rmssO.O 

read(9,*)  n 

write(6,*)  '  n='fn 

read(9t*)  inv 

write(6,*)  *  invs'.inv 

read(9#*)  maxi 

write(6,*)  '  maxi=' ,maxi 

read(9,*)  iwrite 

write(6,*)  '  iwrite: ' , iwrite 


do  100  isl.n 

if(  inv  .eq.  0  )  read(1f2000)  tor(i) 
if(  inv  .eq.  0  )  x(i)=cmplx(  tor(i),0.0  ) 
if (  inv  .eq.  1  )  read(4,2000)  x(i) 
if(  iwrite  .eq.  1  .and.  inv  .eq.  0  )  write(4,*)  x(i) 
100  continue 


call  fft(  x,n,inv  ) 

do  200  1*1, maxi 
mag*cabs(  x(i)  ) 
phase=atan2(  yy,xx  ) 

avemagsamaxl (  cabs(x(1 )) ,cabs(x(2)) ,cabs(x(3)) ,cabs(x(4))  ) 

rmsarms+mag*#2 

write(3,*)  mag/ avemag *100.0 


vv: 


4 


4 


i 


s  * 

V' 


1 


■ 


200 


2000 


continue 

wrlte(6,*)  "averaag=",avemag 

write(6,*)  "rmss",sqrt(  rms*2.0  ) 

write(6,*)  "rma  normalized?", sqrt(  rms/2.0  )/avemag 

format(  elO.M  ) 
call  sclose(O.O) 
stop 
end 


ia  wh  a  w  wsvwswto**  -.'  =? 

hi 

•  .  i 

■>*{ 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


'"  "•V 


c 

c 


subroutine  close 


c 

c 


c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  sclose(dummy) 
closed) 
close(2) 
close (3) 

closed) 
close (7) 
dosed) 
close (9) 
return 
end 
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c**The  following  fft  subroutine  is  taken  almost  verbatim 
c**from  Ahmed  N.,  Rao  K.R.,  Orthogonal  Transforms  for 
^•Digital  Signal  Processing,  Springer-Verlag,  1975.,  pp79 
c 
c 
c 

subroutine  fft(x,n,inv) 
c 

c  this  program  Implements  the  fft  algorithm  to  compute  the  discrete 
c  fourier  coefficients  of  a  data  sequence  of  n  points 
c 

c  calling  sequence  from  the  main  program: 
c  call  fft(x,n,inv) 
c  n:  number  of  data  points 

c  x:  complex  array  containing  the  data  sequence.  In  the  end  dft 

c  coeffs.  are  returned  in  the  array.  Main  program  should 

c  declare  it  as —  complex  x(1024) 

c  inv:  flag  for  inverse 

c  inv=0  for  forward  transform 

c  invsl  for  inverse  transform 

c 

complex  x(1030),w,t,cmplx 
c 

c  calculate  the  #  of  iterations  (log.  n  to  the  base  2) 

c 

itersO 

iremsn 

10  iremslrem/2 

if  (irem.eq.O)  go  to  20 
iter=iter+1 
go  to  10 
20  continue 
slgns-1 

if  (lnv.eq.1)  sign=1 
nxp2sn 

do  50  it si, iter 
c 

computation  for  each  iteration 
nxp:number  of  points  in  a  partition 
c  nxp2:nxp/2 
nxp«nxp2 
nxp2*nxp/2 

wpwr*3. 14 1592/float (nxp2) 
do  40  m«1,nxp2 


calculate  the  multipier 

arg=float(m-1 )*wpwr 

w=cmplx (cos(arg) , sign*sin(arg) ) 

do  40  mxpsnxp,n,nxp 

computation  for  each  partition 

Jl =mxp-nxp+m 
J2=j1+nxp2 
t=x(J1)-x(J2) 
x( jl )sx( jl )+x( J2) 

40  x(j2)st»w 
50  continue 

unscramble  the  bit  reversed  dft 

n2=n/2 
nl =n-1 

J=1 

do  65  isl.nl 
lf(i.ge.j)  go  to  55 
t=x(j) 
x(J)=x(i) 
x(i)=t 
55  k=n2 

60  if(k.ge.j)  go  to  65 
j=j-k 
ksk/2 

go  to  60 
65  j= j+k 

if  (lnv.eq.1)  go  to  75 
do  70  is1,n 
70  x(i)sx(i)/float(n) 

75  continue 
return 
end 


coef fts , 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccoccccccccc 
c  c 

c  subroutine  open  c 

c  c 

ccceccccccccccocccccccccccccccccocccccccccccccccccccccccoccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccoccccccccc 
subroutine  sopen(dummy) 
open(units1,files"torquew) 
rewind ( 1 ) 

open (units2, f ilesndataO/phasett ) 
rewind (2) 

open(units3,fiie=ndataO/Percentn) 
rewind (3) 

open(units4f files"file1") 
rewind (4) 

open(units7(fiiesNdataO/misc2M) 
rewind (7) 

open(units8,filesNdataO/misc3n) 
rewind (8) 

open(units9 , files" input " ) 

rewind (9) 

return 

end 


APPENDIX  D 


PERFORMANCE  PROGRAM 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccocccccccccccccccccccccccc 


c  c 
c  program  to  calculate  developed  power-  c 
c  c 
c  speed  points  for  brushless  dc  motor  c 
c  c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 
c 

c  rpm 

c  wm 

c  p 

c  t 

c  kO 

c  h 

c  w 

c  ws 

c  alpha 

c  ra 

c  xla 

c  rO 

c  xlo 

c  xk 

c  xkg 

c  vd 

c  x 

c  tp 

c  ang 

c  xlamda 

c  m 

c  pd 

c  td 

c  pdl 

c  pd2 

e  pd3 

c  e(1 ) 

c  e<2) 

c  e(3) 

c 
c 

dimension  f(6),v(6),x(6),xx(6),a(6),r(6),e(6)talphs(2) 
common  rr,rd,rO,ra,r(6),pitxla,xlOte(6) ,vm,vdc, ivdon,imodon 
real  iav.irms 
dummy =0.0 
call  popen(dummy) 
pi*3. 1415926536 
eonvapi/180.0 

read(9t*)  istop 


-  speed  of  dc  motor 

-  angular  speed  of  dc  motor 

-  number  of  poles 

-  time 

-  number  of  intervals  between  print  outs 

-  interval  of  integration 

-  electrical  angular  frequency  of  pmm 

-  electrical  angular  frequency  of  pmg 

-  commutation  advance  angle  (degrees) 

-  stator  resistance 

-  stator  leakage  Inductance 

-  choke  coil  resistance 

-  choke  coil  Inductance 

-  motor  constant 

-  generator  constant 

-  inverter  input  voltage 

-  array  of  stator  phase  currents 

-  period 

-  electrical  angular  displacement 

-  electrical  angular  displacement  (pmg) 

-  number  of  phases 

-  total  developed  power 

-  total  developed  torque 

-  phase  1  developed  power 

-  phase  2  developed  power 

-  phase  3  developed  power 

-  excitation  voltage  of  phase  1 

-  excitation  voltage  of  phase  2 

-  excitation  voltage  of  phase  3 
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write<6,»)  *  istop* istop 

read (9,*)  alpha 

write<6,»)  ’  alphas’, alpha 

read(9,*)  (  a(i),isl,5  ) 

write(6,«)  ’  as»,(  a(i),isi,5  ) 

read(9,*)  h 

write(6,#)  ’  h=’,h 

read(9,*)  rO 

write(6,»)  *  rOs*,rO 

read(9»*)  xlO 

write(6,»)  ’  xlOs’.xlO 

read(9,*)  ra 

write(6,»)  ’  ras’.ra 

read(9,*)  xla 

write(6,»)  *  xlas'.xla 

read(9,*)  rd 

write(6, *)  ’  rds’.rd 

read(9,*)  rpmO 

write(6,*)  ’  rpmO=’,rpmO 

read (9,*)  xk 

write(6,#)  ’  xk=’,xk 

read(9,*)  xkg 

write (6,*)  '  xkg=' ,xkg 

read(9,*)  poax 

write(6,»)  '  pmaxs'.pnax 

read(9,*)  pmin 

write(6,*)  •  pmlns'tpmin 

read(9,*)  vdc 

write (6,*)  *  vdc=',vdc 

read (9,*)  points 

write(6,»)  •  pointss*, points 

read(9»*)  ivdon 

write (6,*)  ’  ivdons ivdon 

read(9,*)  imodon 

write(6,*)  '  imodon s’ , imodon 

read (9,*)  ichon 

wrlte(6,»)  '  ichons', ichon 

read(9,*)  lwgon 

write (6,*)  '  iwgons ’ , lwgon 

read(9,»)  (  x(i),isi,3  ) 

write(6,»)  ’  xs’,(  x(i),isi,3  ) 

read (9,*)  iarea 

write(6,*)  *  iareas iarea 

read(9,*)  xwg 

write (6,*)  ’  xwgs ’ , xwg 

read(9,*)  ifull 

write(6,»)  '  ifulls’, ifull 

do  4  Js1,5 

if(  Imodon  .eq.  1  )  a(j)s  a( J)*pi/180.0 
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If (  imodon  .eq.  0  )  a(j)=pi/2.0 
if(  imodon  .eq.  0  )  a(1 )s30.00*pi/180.0 
4  continue 

m=3 

pfw= (rpm0/45000 . 0 )«*2»330 . 0 

wm=rpm0*pi/30.0 

if(  iwgon  .eq.  1  )  wg=7950*pi 

if(  iwgon  .eq.  0  )  wg=xwg*wm 

wmsabs(wm) 

p=4.0 

t=0.0 

w=p/2.0*wm 

wa=2.0*wg 

vm= sqrt C 3 . 0 ) * xkg*wg 
write(6,*)  *  vats', vm 
alpha=alpha*conv 

areaa0.5+0.5*(  pi/3.0-a(1)+a(2)-a(3)+a(4)-a(5)  )/(  pi/3.0  ) 

if C  imodon  .ne.  0  .and.  iarea  .eq.  1  )  vdcsvdc/area 

if(  ifull  .eq.  1  )  alphs(1 )=acoa(  (vdc#sqrt(2.0))/(1.0»vm#1.35)  ) 

if(  iftill  .eq.  1  )  alpha(2)=acoa(  (vde#aqrt(2.0))/(2.0*vm<,1.35)  ) 

if(  ifull  .eq.  0  )  alphsd )»acos(  (vdc#sqrt(2.0))/(0.5*vm»1.35)  ) 

if(  ifull  .eq.  0  )  alpha(2)=acoa(  (vdc#sqrt(2.0))/(1.0»vmB1.35)  ) 

write(6,»)  *  alphad  )='  ,alpha(1  )/eonv 

write(6,*)  *  alphs(2)s',alpha(2)/conv 

write(6,*)  ’  areas*, area 

write(6,»)  '  fs» ,w/(2.0#pi) 

write(6,»)  *  fas' ,ws/(2.0*pi) 

write(6,*)  '  mag  of  els' ,xk*w/2.0 

xxd  )=0.0 

xx(2)a0.0 

xx(3)=0.0 

pdsO.O 

pinsO.O 

iavsO.O 

irmaaO.O 

elavesO.O 

vdavesO.O 

e IrmaaO.O 

sumaO.O 

tp*2.0*pi/w 

tminapmln*tp 

tmax*pmax*tp 

tmxmnatmax-tmin 

scales30.0/pointa 

kOaint(  tmxmn/tp*acale/aba(rpmO)/h  ) 
kountsO 
100  index* 1 

rewind (9) 
read(9#*)  iatop 
if(  iatop  .eq.  1  )  go  to  305 
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kountskount+1 
200  ang=w*t 

ipr in 1 sint ( ang/2 . 0/pi ) 
angaang-floatUprinl  )*2.0*pi 
e(1 )=xk*w*2.0/p*sin(ang-alpha) 
e(2)sxk*w*2.0/p*sin(ang-alpha-2.0/3.0*pi) 
e ( 3 )=xk*w*2 . 0/p*sin ( ang-al pha+2 . 0/3 . 0#pi ) 
wstawa*t0 

if (  wst  .ge.  2.0*pi  )  tOsO.O 
wat=wa*t0 

if (  ifull  .eq.  0  )  call  half(  wat, xlamda  ) 
if(  ifull  .eq.  1  )  call  full(  wst , xlamda  ) 
call  fv(  a , ang,v,vd, alpha, xlamda  ) 
call  ff(  index, t, f,x,xi , ang, del, v  ) 
if(  ichon  .eq.  1  )  call  checkC  m,x,f  ) 
call  derk(  m,x,f,t,tO,h, index  ) 
lf(  index  .eq.  2  )  go  to  200 
if(  t  .It.  tmin  )  go  to  100 
f ( 1 )* (  x(1)-xx(1)  )/h 
f (2)=(  x(2)-xx(2)  )/h 
f(3)=<  x(3)-xx(3)  )/h 
do  86  kks 1 , 3 
xx(kk)=x(kk) 

86  continue 

call  ffi(  fi,f ,ang  ) 

call  fvap(  wat,vap,alphaa  ) 

v 1 p* vd-f i *xl0/2 . 0-x i*r0-x ( 1 ) • ( r ( 1 ) -ra ) 

vacravap-vlp 

if(  ang  .It.  3.0/3.0*pl  .and.  wat  .It.  2.0/3.0*pi  ) 
c  vcarax(1)*(r(1)-ra) 

lf(  ang  .ge.  3.0/3.0*pi  .and.  wat  .ge.  3.0/3.0*pl  .and. 
c  wat  .It.  4.0/3.0*pi  )  vacrax(1)#(r(1)-ra) 

85  if(  t  . ge.tmax  )  go  to  300 

lf(  kount  .It.  kO  )  go  to  100 

300  sums sum* 1 . o 


pd1ax(1)*e(1) 

pd2ax(2)*e(2) 

pd3*x(3)*e(3) 

tda(  pd1+pd2+pd3  )/wo 

pd=pd+(  pd1+pd2+pd3  ) 

pinapin+(  v(1 )»x(1 )+v(2)*x(2)+v(3)#x(3)  ) 

irmaairna+x(1)*B2 

iavsiav+abs(  x(1)  ) 

e 1 rna*e 1 raa+e ( 1 ) #*2 

elaveaelave-fabaC  e(1)  ) 

vdave*vdave-«>vd 

hsin6ahsin6+td *ain ( 6 . 0*ang ) 

hcoa6ahcos6+td *coa ( 6 . 0*ang ) 

hain12ahain12+td#ain(12.0fang) 

hcoa 1 2 ahcoa 1 2+td *coa(12.0Bang) 


hsin18=hsin18+td*sin(18.0*ang) 

hcos 1 8=hcoa 1 8+td  *cos ( 1 8 . 0*ang ) 

hsin24shaln24+td*sin(24.0»ang) 

hcos24=hcos24+td*coa ( 24 . 0*ang ) 

txat-tmxmn 

txxat*w* 180. 0/pi 

write ( 1 , 2000 )  x(1) 

wrlte(2,3000)  txx 

write (3, 2000)  td 

wrlte(4,2000)  v<1) 

write (7 ,2000)  vd 

write (8, 2000)  x(1 )+x(2)+x(3) 

kountsO 

if(  t  .It.  tmax  )  go  to  100 
305  pdspd/suo 

pinapin/aum 
Irma* Irma /sum 
lrmsssqrt(  Irma  ) 
lavsiav/sum 
e 1 rmsae 1 rms/sum 
e1rma=sqrt(  elrma  ) 
elavese lave/sum 
vdavesvdave/sum 
hsln6=haln6/sum 
hcos6=hcos6/sum 
h6=hsin6*#2+hcos 6*»2 
h6ssqrt(  h6  )/6.0 
haln12=hain 12/sum 
hcos 1 2=hcos 1 2/ sum 
h  1 2  *hs  1  n  1 2M2+hcos  1 2  **2 
h12aaqrt(  h12  )/12.0 
hslnl8ahsin 18/sum 
hooa 1 8=hcos 1 8/sum 
h  I8ahsln  l8**2+hcos  1 8”2 
h18=aqrt(  h18  )/18.0 
hsln24ahsin24/sum 
hcos24shcos24/aum 
h24ahsin24*»2+hcos24*»2 
h24asqrt(  h24  )/24.0 
600  tsava ( pd-pfw ) /wm 

write (6,*)  '  tsava ',taav 
effa(pd-pfw)/pln*100.0 
lf<  eff  .It.  100.0  )  go  to  710 
effapin/(pd+pfw)B 100.0 
710  write(6,*)  '  pina'.pin 

write(6,»)  '  effa'.eff 
wrlte(6,*)  '  lava', lav 
wrlte<6,»)  '  lrmsa',irms 
wrlte(6,*)  *  elavea' , slave 
wrlte(6,f)  '  vdavea' ,vdave 
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V-  rf 


*  Ajt-  %  X  VI 4  >  ^LriA  |V  V 


L4tU  4L#U  *.I 


1 


write(6,»)  •  h6='#h6 
write(6,»)  •  h12s',h12 
write(6,»)  »  h18=* ,h18 
write(6,»)  •  h24s»,h24 
formate  5(e12.4,5x)  ) 
formate  e10.4  ) 
formate  f8.3  ) 
call  pcloseedummy) 
stop 
end 


-VI 


t>| 


I 


cccccccocccocecccccccecccccccccccccccecccccccccccccccccccccccccccccccccc 
ccecccccccceccccccccccccccccccccccccccccccccccccccccccccccccccccccceccec 
c  c 

c  subroutine  check  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  check(  m,x,f  ) 
dimension  x(6),f(6) 

common  rr,rd,r0,ra,r(6),pi,xla,xl0, e(6),vmtvdc,ivdon,imodon 
do  10  i=1 ,m-1 

slgnf=abs(  f(i)  )/f(i) 
signx=abs(  x(l)  )/x(l) 

lf(  abs(  f(i)  )  .gt.  1.0e08  )  f (i)ssignf»1 .0e08 
lf(  abs(  x(l)  )  .gt.  1.0e06  )  x(i)=signx*1.0e06 
10  continue 

f(3)*-(  f(1)+f(2)  ) 

return 

end 


«L  1L 


cceccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccocccccccccccccccecccccccecccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  pclose  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccccccc 
cecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  pclose (dummy) 
close ( 1 ) 
close (2) 
close (3) 

close(4) 
close (7) 
dose(8) 
close (9) 
return 
end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  derk  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  derk  is  a  fourth-order,  fixed  increment  runge-kutta 

c  integration  routine 

c 

cm-  number  of  simultaneous  differential  equations 
c  x  -  array  of  dependent  variables 

c  f  array  of  derivatives  of  independent  varibles 

c  h  -  increment 

c  t  time 

c  index  -  indicator 

c  1.  if  exit  routine  with  index=1,  solution  found  atst+h 

c  2.  if  exit  routine  with  indexa2,  go  back  to  re-evaulate 

c 

subroutine  derk(  m,x,f ,t,tO,h, index  ) 
dimension  x(6)(f(6),q(400) 
lf(  index  .eq.  2  .)  go  to  19 

18  kxxsO 
index=2 

do  35  i=1 ,m 
j=i+300 

35  qCj)sx(i) 

19  kxx=kxx+1 

go  to  (1,2,3,4),kxx 

1  do  5  i=1 ,m 
q(i)af(i)*h 

5  x(i)sx(i)+q(i)/2.00 
tat+h/2.00 
t0at0+h/2.00 
return 

2  do  6  ia1,m 
JaUlOO 
kai+300 
q(J)af(i)»h 

6  x(i)aq(k)+q(J)/2.00 
return 

3  do  7  ia1,m 
Jai+200 
kai+300 
q(J)af(i)»h 

7  x(i)sq(k)+q(J) 
t*t+h/2. 00 
t0at0t*/2.00 
return 
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yf 


S 
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& 


■«r. 

«,V] 


do  8  i=1 ,m 
jsi+100 
kai+200 
1=1+300 

x(i)=q(l)+(  (  q(i)+2.00*q(k)+2.00*q(k)+f(i)*h  )/6.00  ) 

index =1 

return 

end 


cccccccccccccccccccccccccccocccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  c 

c  subroutine  f  c 

c  c 

cccccccccccccccccccccccccccccccceccccceccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccoccccccccccccccccccccccccccccccccccccccccccccccccc 


subroutine  ff(  index,t,f ,x,xi,ang,del,v  ) 
dimension  f (6) ,v(6) ,x(6) 

common  rr,rd,r0,ra,r(6) ,pi,xla,xl0,e(6) ,vm, vdc , ivdon , imodon 
rrsrd+ra 


if(  index  .eq.  1  )  r(1)=rr 
if(  index  .eq.  1  )  r(2)srr 
if(  index  .eq.  1  )  r(3)=rr 


if  ( 

ang 

.ge. 

0. 

0 

.and. 

ang 

.It. 

1.0/3.0*pi 

) 

goto 

10 

if  ( 

ang 

.ge. 

1. 

0/3.0»pi 

.and. 

ang 

.It. 

2.0/3.0*pi 

) 

goto 

20 

if  ( 

ang 

.ge. 

2. 

0/3.0«pi 

.and. 

ang 

.It. 

3.0/3.0»pi 

) 

goto 

30 

if( 

ang 

.ge. 

3. 

0/3.0»pi 

.and. 

ang 

.It. 

4.0/3.0*pi 

) 

goto 

40 

if( 

ang 

.ge. 

4. 

0/3.0«pi 

.and. 

ang 

.It. 

5.0/3.0»pi 

) 

goto 

50 

if( 

ang 

.ge. 

5. 

0/3.0»pi 

.and. 

ang 

.It. 

6.0/3.0»pi 

) 

goto 

60 

10  n1=1 

n2=3 
n3=2 

if(  icycle  .eq.  1  )  lflagsO 
i tests 10 
icyclesO 
go  to  100 

20  nisi 

n2s2 

n3*3 

itests20 
go  to  200 
30  n1s2 

n2si 

n3*3 
itests30 
go  to  100 
40  n1s2 

n2s3 

n3*1 

itests40 
go  to  200 
n1s3 
n2s2 
n3*1 

itestsSO 
go  to  100 


50 


n2=1 
n3=2 
iteat=60 
lcyolesl 
go  to  200 

100  lf(  x(n2)  .It.  -0.001  .and.  index  .eq.  1  )  iflag=itest 

if(  x(nl)  .It.  -0.001  .and.  index  .eq.  1  )  r(n1 )srr+1000.0*abs(  x(nl)  ) 

if(  x(n3)  .gt.  0.001  .and.  index  .eq.  1  )  r(n3)=rr+1000.0*abs(  x(n3)  ) 

if(  iflag  .eq.  itest  .and.  index  .eq.  1  )  r(n2)=rr+1000.0*abs(  x(n2)  ) 

do  150  is1,3 

if (  r(i)  .gt.  1000.0  )  r(i)=1000.0 
150  continue 

xlsxla 

f(n1)*(  -x(n1 )*r(n1 )+v(n1 )-e(n1 )  )/xl 
f(n2)=(  -x(n2)*r(n2)+v(n2)-e(n2)  )/xl 
f(n3)=(  -X(n3)*r(n3)+v(n3)-e(n3)  )/xl 
return 

200  if(  x(n2)  .gt.  0.001  .and.  index  .eq.  1  )  iflag= itest 

if(  x(n3)  .gt.  0.001  .and.  index  .eq.  1  )  r(n3)arr+1000.0»abs(  x(n3)  ) 

if(  x(nl)  .It.  -0.001  .and.  index  .eq.  1  )  r(n1)=rr-*-1000.0«abs(  xCnl)  ) 

if(  iflag  .eq.  itest  .and.  index  .eq.  1  )  r(n2)=rr+1000.0*abs(  x(n2)  ) 

do  250  1=1,3 

if(  r(i)  .gt.  1000.0  )  r (i)=1000.0 
250  continue 

xlsxla 

f(n1)=(  -x(n1)*r(n1)+v(n1)-e(n1)  )/xl 
f(n2)=(  -x(n2)*r(n2)+v(n2)-e(n2)  )/xl 
f(n3)*(  -x(n3)*r(n3)+v(n3)-e(n3)  )/xl 


return 

end 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccecceccccccecccccccccccccccecccccccccccccccccccccccccccccccecccccccc 
c  c 

c  subroutine  fi  c 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccceccccc 
ccccccccccoeececcccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  ffi(  fi,f,ang  ) 
dimension  f(6) 

common  rr,rd,r0,ra,r<6) ,pi,xla,xlO, e(6) ,vm,vdc,ivdon,imodon 


ge. 

0, 

,0/3.0«pi 

.and. 

ang 

.It. 

1.0/3.0«pi 

) 

fi=-f (2) 

ge. 

1. 

.0/3.0*pi 

.and. 

ang 

.It. 

2.0/3.0*pi 

) 

fi= 

f(1) 

ge. 

2. 

,0/3.0*pi 

.and. 

ang 

.It. 

3.0/3.0*pi 

) 

fi=-f (3) 

ge. 

3. 

,0/3.0»pi 

.and. 

ang 

.It. 

4.0/3.0*pi 

) 

fis 

f(2) 

ge. 

4. 

,0/3.0«pi 

.and. 

ang 

.It. 

5.0/3.0*pi 

) 

fi=-f(1) 

ge. 

5. 

,0/3.0«pi 

.and. 

ang 

.It. 

6.0/3.0*pi 

) 

fis 

f  (3) 

return 

end 


cccccccccccccccccocccccccccccccecccceccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccecccccccccccccccecccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  full-bridge  c 

c  c 

cccccccceccccccccccccecccccccoccccecoccccccccccccccocccccccccccccccccccc 
cecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccccccccc 

subroutine  full(  wst.xlamda  ) 

common  rr(rd,rO,ra,r(6),pi(xla,xlOte(6),vm,vdc,ivdon,imodon 

if(  wst  .ge.  0.0  .and.  wst  .It.  1.0*pi/3.0  )  xlamdas 

c  wst-0.0 

if(  wst  .ge.  1.0*pi/3.0  .and.  wst  .It.  2.0*pi/3.0  )  xlamdas 
c  wst-1.0*pi/3.0 

if(  wst  .ge.  2.0*pi/3.0  .and.  wst  .It.  3.0*pi/3.0  )  xlamdas 
c  wst-2.0*pi/3.0 

if(  wst  .ge.  3.0*pi/3.0  .and.  wst  .It.  4.0*pi/3.0  )  xlamdas 
c  wst-3.0*pi/3.0 

lf(  wst  .ge.  4.0*pi/3.0  .and.  wst  .It.  5.0*pi/3.0  )  xlamdas 
c  wst-4.0*pi/3.0 

ifC  wst  .ge.  5.0vpi/3.0  .and.  wst  .It.  6.0*pi/3.0  )  xlamdas 
c  wst-5.0*pi/3.0 

return 

end 
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cceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


subroutine  half-bridge 


cccccccccccccccccecccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccecceccccccccccccccccccccccccccccccccccccccc 


subroutine  half(  wst.xlarada  ) 

common  rrfrd,r0,ra,r(6),pitxla,xl0,e(6) ,vm, vdc.ivdon, iraodon 
if(  wst  .ge.  0.0  .and.  wst  .It.  2.0*pi/3.0  )  xlamda= 

c  wst-0.0 

if (  wst  .ge.  2.0*pi/3.0  .and.  wst  .It.  4.0*pi/3.0  )  xlamda= 
c  wst-2.0*pi/3.0 

ifC  wst  .ge.  4.0*pi/3.0  .and.  wst  .It.  6.0»pi/3.0  )  xlamda= 
c  wst-4.0*pi/3.0 
return 
end 
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cccccccceececccccccecccceccccccccccccccccccccccccccccccccecccecccccccccc 
ccccccccececcecccccccccccccccccecccccccceecccccccccccccccccccccccceccccc 
c  c 

c  subroutine  popen  e 

c  c 

eecccccceeeccccccccccecccccecccccccccccccccccccccccccocccccccceccceccccc 
ccccccccccccccccccoeccceoocoocecceccccccecceccccccceooeccccccccccccccccc 
subroutine  popen (dunmy) 
open (unit s1,files"dataO/I1") 
rewlnd(l) 

open (units2, f ilesNdataO/Tiae" ) 
rewind (2) 

open(units3ffiles"dataO/TorqueH) 
rewind (3) 

open(units4yfile3,,dataO/Van,,) 
rewind (4) 

open (unit =7 , file* "dataO/Vd" ) 
rewind (7) 

open (unit s8tflle3”dataO/In") 
rewind (8) 

open (unit =9 ,  files" input " ) 

rewind (9) 

return 

end 


cceccccccccccceccceccoccccccccecccccccceeeeecccccccccccccccccccccccccccc 
cccccccoceeccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  subroutine  vd  c 

o  c 

cceeccceocccocccccceeocceccccccoceeeeececccccccccccccccccccecccccccccccc 
ccccccccccececcccecccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  fv(  a,angfv,vd,alphs,xlamda  ) 
dimension  a(6),v(6),alphs(2) 

common  rr,rd,r0,ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon,imodon 
beta=xlamda+alphs ( 1 ) 

if(  ivdon  .eq.  1  .and.  ifull  .eq.  1  >  vd=vm*sin(  beta+pi/3.0  ) 

if (  ivdon  .eq.  1  .and.  ifull  .eq.  0  )  vd=vm*sin(  beta+pi/6.0  ) 

if(  ivdon  .eq.  0  )  vdsvdc 

shift=30. 0/180. 0*pi 

theta=shift+ang 

call  svxn(  vx, theta, a  ) 

vansvx 

thetasshif t+ang-2. 0/3 . 0*pi 
call  svxn(  vx, theta, a  ) 
vbn«vx 

thetaxshift+ang+2.0/3.0*pi 
oall  svxn(  vx, theta, a  ) 
vcnxvx 

v(1)«vd*van 

v(2)svd»vbn 

v(3)*vd#vcn 

return 

end 


cccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccccccc 
ecceecccceeeecceceecceccccocccccccccceccccccccccccccoccccccececccccccccc 
c  c 
e  subroutine  fvap  c 
c  c 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cceceococoooeocccooccooccoeeeoeeceocooccooococececccccecccccoeccceccccec 


subroutine  fvap(  wst,vap,alps  ) 

common  rrtrd,rOtratr(6),pl(xlalxlO,e(6)lvm,vdc,lvdon,lmodon 
vaasO.O 

vabs  vm*sln(  wst+1.0/3.0*pl+alpa  ) 
vac=-vm*sin(  wst+5.0/3.0*pl+alps  ) 

if(  wst  .ge.  0.0/3*0*pl  .and.  wst  .It.  1.0/3.0*pi  )  vap*  vab 

lf(  wst  .ge.  1.0/3.0*pl  .and.  wst  .It.  2.0/3. 0*pl  )  vap*  vac 

lf(  wst  .ge.  2.0/3.0*pl  .and.  wst  .It.  3.0/3.0*pi  )  vap*  vac 

lf(  wst  .ge.  3.0/3.0*pl  .and.  wst  .It.  4.0/3.0«pl  )  vap*  vaa 

lf(  wst  .ge.  4.0/3.0*pl  .and.  wst  .It.  5.0/3.0*pl  )  vap*  vaa 

lf(  wst  .ge.  5.0/3.0*pi  .and.  wst  .It.  6.0/3.0*pl  )  vap*  vab 

return 
end 
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ccecccccccccccccececcccecccccececccecccccccccccccccccccccccccccccccccccc 
ccccccecccccccccccccccccccccccccccccccccccccccccccccccccceccccccccccccec 
c  c 

c  subroutine  vxn  c 

c  c 

cccececcccccccccccceccccccccccccceecccccccccccccceeccccccccccccccccccccc 
ccccccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  svxnC  vx,th,a  ) 
dimension  a(6) 

common  rr,rd,rO,ra,r(6),pi,xla,xlO,e(6) ,vm, vdc,ivdon,imodon 

vxsO.O 

xa=0.0 

xbspi 

xc=2.0*pi 

a(6)=pi/2.0 

th=th+2.0*pi 

th=th-float(  int(th/xc)  )*xc 

if(  th  .ge.  xa+a(1)  .and.  th  .It.  xa+a(2)  )  vx=1.0 

if(  th  .ge.  xa+a(3)  .and.  th  .It.  xa+a(4)  )  vx=1.0 

if(  th  .ge.  xa+a(5)  .and.  th  .It.  xa+a(6)  )  vx=1.0 

if(  th  .ge.  xb-a(6)  .and.  th  .It.  xb-a(5)  )  vx=1.0 

if(  th  .ge.  xb-a(4)  .and.  th  .It.  xb-a(3)  )  vxsl.0 

if(  th  .ge.  xb-a(2)  .and.  th  .It.  xb-a(1)  )  vx=1.0 

lf(  th  .ge.  xb«>a(1)  .and.  th  .It.  xb+a(2)  )  vx=-1.0 

lf(  th  .ge.  xb-t-a(3)  .and.  th  .It.  xb-<>a(4)  )  vx*-1.0 

if(  th  .ge.  xb-»-a(5)  .and.  th  .It.  xh«>a(6)  )  vx*-1.0 

if(  th  .ge.  xc>a(6)  .and.  th  .It.  xc>a(5)  )  vx*-1.0 

lf(  th  .ge.  xc-a(4)  .and.  th  .It.  xc-a(3)  )  vx*-1.0 

if(  th  .ge.  xc-a(2)  .and.  th  .It.  xc-a(1)  )  vx*-1.0 

return 
end 


